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FOREWORD 

This  semiannual  technical  report  entitled,  "Applica- 
tion of  Advanced  Methods  for  Identification  and  Detection  of 
Nuclear  Explosions  from  the  Asian  Continent,"  is  submitted 
by  Systems,  Science  and  Software  (S3)  to  the  Advanced  Research 
Projects  Agency  and  to  the  Air  Force  Office  of  Scientific 
Research  (AFOSR) . This  report  presents  the  results  of  a con- 
tinuing effort  to  obtain  an  optimum  multi-discriminant/ 
detection  procedure  foi  earthquakes  and  explosions  occurring 
within  the  Asian  Continent.  The  work  is  being  performed 
under  Contract  Number  F44620-74-C-0063.  Mr.  William  J.  Best 
is  the  AFOSR  technical  contracting  officer. 

Dr.  J.  Theodore  Cherry  is  the  S3  project  manager.  Drs. 
Thomas  C.  Bache  and  Joseph  F.  Masso  are  responsible  for  the 
development  and  application  of  the  seismic  ground  motion  pre- 
diction work.  Dr.  John  M.  Savino  and  Mr.  Kenneth  G.  Hamilton 
are  responsible  for  the  analysis  of  the  seismic  data  against 
which  all  theoretical  developments  must  eventually  be  tested. 
Acting  as  consultants  on  the  project  are  Professors  Charles 
B.  Archambeau  of  the  University  of  Colorado,  David  G. 

Harkrider  of  the  California  Institute  of  Technology  and 
Donald  V.  Helmberger  of  the  California  Institute  of  Technology. 

The  authors  wish  to  extend  their  sincere  appreciation 
to  Ms.  Bernadine  Ludwig  and  Ms.  Darlene  Roddy  for  the  many 
hours  spent  on  the  preparation  of  this  report. 
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I.  INTRODUCTION 


The  fundamental  problem  toward  which  our  rose,' 


""ivii  uui  n'bearcn  nas 

een  directed  is  seismic  discrimination.  That  is,  from  the 
background  of  an  earthquake  and  microseismic  signals,  one 

wishes  to  identify  those  teleseismic  events  which  are  caused 
by  nuclear  explosions. 


Our  approach  to  the  problem's  solution  has  been  the 
development  of  a deterministic  methodology  for  predicting 
teleseismic  ground  motion  from  both  earthquakes  and  nuclear 
explosions.  This  predictive  capability  provides  a theoretical 
basis  against  which  existing  discriminants  may  be  tested  and 
new  discriminants  designed. 


Figure  1.1  is  a schematic  of  the  technique  used  to 
predict  the  body  wave  seismogram  generated  by  a nuclear  explo- 
sion. This  model  has  been  separated  into  five  parts,  which 
are  treated  independently.  These  are: 


1.  Calculation  of  the  equivalent  elastic  source 
generated  by  the  explosion. 

2.  Modification  of  the  equivalent  source  by  re- 
verberations in  the  near  source  crustal  struc- 
ture . 


3.  Propagation  of  the  rays,  emerging  from  the  base 
of  the  crust,  through  the  upper  mantle  to  the 

base  of  the  crust  in  the  vicinity  of  the  re- 
ceiver . 

4.  Pulse  modifications  by  crustal  reverberations 
at  the  receiver. 

5.  Modification  of  the  computed  ground  motion  by 
the  characteristics  of  the  seismometer. 

In  a similar  fashion,  given  the  equivalent  elastic 
source  from  the  explosion,  teleseismic  ground  motion  from 
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Rayleigh  and  Love  surface  waves  are  computed  using  a computa- 
tional technique  based  on  the  Thomson-Haskell  matrix  method 
(Harkrider  (1964)). 

Both  body  wave  and  surface  wave  teleseismic  ground 
motion  from  earthquakes  may  be  predicted  if  the  earthquake 
equivalent  source  replaces  the  explosion  equivalent  source. 

Therefore,  our  research  may  be  subdivided  into  the 
following  four  categories: 

1.  The  prediction  of  the  equivalent  elastic  source 
from  nuclear  explosions. 

2.  The  prediction  of  the  equivalent  elastic  source 
from  earthquakes. 

3.  Propagation  of  those  sources  to  teleseismic  dis- 
tances including  both  body  waves  and  surface 
waves . 

4.  Evaluation  of  existing  discriminants  and  identi- 
fication of  improved  discriminants  based  on 
these  predictions. 

For  the  most  part  the  predictive  capability  required 
for  the  analysis  (Items  1,  2 and  3)  is  now  available.  These 
theoretical  models  are  being  used  to  predict  the  equivalent 
elastic  sources  from  earthquakes  and  explosions  and  to 
propagate  these  sources  to  teleseismic  distances.  Further- 
more, based  on  expected  differences  between  earthquake  and 
explosion  sources,  a new  discriminant  is  being  evaluated. 

Section  II  of  this  report  presents  our  current  under- 
standing of  the  equivalent  elastic  source  generated  by  a 
nuclear  explosion.  These  source  functions  represent  the 
coupling  of  the  explosion  into  the  elastic  regime,  i.e., 
seismic  coupling.  The  effects  of  elastic  and  inelastic  near 
source  material  properties  and  changes  in  explosion  yield  are 
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identified  and  interpreted  in  terms  of  changes 
magnitude . 


in  teleseism ic 


In 

computed 

functions 

and  on  m, 
o 

explosion 
tectonic 
the  crust 


Section  III,  theoretical  body  wave  seismograms  are 
at  teleseismic  distances  using  specific  source 
from  Section  II.  The  effect  on  these  seismograms 
is  presented  as  a function  of  depth  of  burial, 
yield,  the  inelastic  properties  of  the  upper  mantl 
stress  release  and  near  source  reverberations  in 


e , 


Section  IV  presents  the  results  of  a variable  frequency 
magnitude  (VFM)  discriminant  operating  on  body  wave  seismo- 
grams which  exploits  theoretically  predicted  spectral  dif- 
ferences in  the  equivalent  source  between  earthquakes  and 
explosions.  Complete  separatim  of  earthquake  and  explosion 
populations  is  achieved  down  to  body  wave  magnitudes,  m, 
of  the  order  4.0  to  4.5. 


V 


In  Section  V a comparison  is  made  between  the  effec- 
tiveness of  the  VFM  discriminant  and  other  proposed  dis- 
criminants (complexity,  spectral  ratio,  higher  moments  of 
frequency  and  mb  - Mg) . All  discriminants  are  applied  to 
the  same  data  base:  the  large  Eurasian  population  recorded 

at  LASA.  It  appears  that  a combined  application  of  the  VFM 
and  modified  moment  of  frequency  (MMF)  discriminants  may  pro- 
vide a very  effective  technique  for  event  identification. 


The  appendices  (A  through  G)  present  theoretical 
justification  for  a number  of  key  issues  presented  in  the 
various  sections  of  the  report.  Sections  IV  and  V should  be 
read  in  sequence.  All  other  sections  and  the  appendices  are 
essentially  self-contained  and  may  read  independently. 
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II.  THE  EQUIVALENT  ELASTIC  SOURCE  FROM  A NUCLEAR  EXPLOSION 

2.1  INTRODUCTION 

The  common  element  required  for  predicting  both  body 
wave  and  surface  wave  teleseismic  ground  motion  from  a nuclear 
explosion  is  the  equivalent  elastic  source  of  the  explosion. 

By  and  large,  seismologists  have  been  forced  to  use 
experimentally  determined  explosion  source  functions  for 
teleseismic  ground  motion  predictions.  These  source  func- 
tions are  limited  in  number,  quality  and  the  rock  environ- 
ment they  represent.  These  limitations  are  due  mainly  to 
the  expense  and  technical  difficulty  involved  in  fielding  an 
experiment  designed  to  obtain  free  field,  time  history  data 
in  the  elastic  regime  for  a nuclear  source. 

Computational  techniques  and  material  models  have  been 
developed  at  S3  which  permit  a prediction  of  the  equivalent 
elastic  source  given  the  yield  of  the  explosive  and  the 
material  properties  of  the  near  source  geologic  environment. 

Since  laboratory  tests  on  rock  samples  provide  the 
essential  data  base  for  a preshot  estimate  of  material  pro- 
perties, it  is  important  to  determine  the  sensitivity  of  the 
equivalent  source  to  those  material  properties  which  are 
subject  to  routine  laboratory  measurement.  In  this  chapter 
we  present  the  results  of  a computational  parameter  study 
which  identifies  the  dependence  of  teleseismic  magnitude  on 
near  source  material  properties. 

2 . 2 THE  REDUCED  DISPLACEMENT  POTENTIAL 

If  u(r,t)  is  the  radial  displacement  from  a spherically 
symmetric  explosion  (see  Appendix  G for  the  general  equivalent 
source  representation)  where  u(r,t)  is  measured  in  the  linear 
elastic  region,  then 
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(2.1) 


where 


Ur)  _ 

p ■'  - the  equivalent 


source  potential 


^Ct)  the  reduced  displacement  potential  (RDP) 
t = t - r/a 

a = P wave  velocity  in  the  source  region 

For  body  waves  we  save  only  the  r'1  component.-  in  the 
equivalent  source  expansion.  Therefore 


u(r,t)  c lill 
ra 


(2.2) 


For  surface  waves  (see  Appendix  D for  a detailed  analysis 
of  Rayleigh  wave  amplitude  scaling)  the  forcing  function  is  an 

equivalent  radial  stress  orr  applied  to  the  interior  of  a spheri 
cal  cavity  where 


rr 


iiil  + 4g2  / iLii  nD\ 

T \ r3  r2a  / 


(2.3) 


where  p is  the  density 

$ = Ai/p  is  the  shear  wave  velocity 
P is  the  shear  modulus. 

For  those  frequencies  of  interest  for  teleseismic 
surface  waves,  Eq.  (2.3)  reduces  to 


rr 


4y  Hi! 


(2.4) 
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Equations  (2.2)  and  (2.4)  isolate  the  elastic  properties 
of  the  source  region  responsible  for  body  wave  and  surface 
wave  amplitude.  Body  wave  amplitude  is  inversely  proportional 
to  compressional  wave  velocity  a,  while  surface  wave  ampli- 
tude is  directly  proportional  to  the  shear  modulus  y. 

If  we  assume  that  over  a limited  frequency  band 


-n 


(2.5) 


i. 


where  V is  the  amplitude  spectrum  of  T and  f is 
frequency,  then  at  a given  yield,  W , body  wave  amplitude 
is  proportional  to 


V 


af 


n 


(2.6) 


where  f is  the  dominant  frequency  of  the  pulse.  If  the 

combination  of  the  earth’s  Q and  the  seismometer  band  pass 

is  very  narrow  with  center  frequency  f then  at  a new 

yield  W 0 

2 


¥ = 


A W 

A 2 


af 


•n  W 


where  f cube  root  scales  into  f , i.e., 
* o 


/W  \1/3 

" M<) 


Therefore 


'v  W 

• A * 

XU  r=  M 2 

2 *n  W 
afe  1 


(l-n/3) 


= Y 


(2.7) 
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Fcr  surface  waves,  Eqs . (2.4)  and  (2.6)  give 


A 

Y 


U 


A 

f (Vi) 


0 


o 


Changing  the  yield  to  W2  and  assuming  a narrow  band 
filter  gives 


A 


Y 

2 


(2.8) 


This  scaling  law  is  identical  to  that  obtained  for  body 
waves  (Eq.  (2.7)).  Both  surface  waves  and  body  waves  scale 
with  yield  according  to^the  l-n/3  power  where  n corres- 
ponds to  the  slope  of  Y in  the  frequency  band  of  interest 

If  the  source  spectrum  is  flat,  then  n = 0 and  tele- 
seismic  amplitudes  scale  with  the  first  power  of  the  yield. 


2.3  CALCULATION  OF  THE  RDP 

If  we  assume  that  the  displacement  is  linear  within  a 
given  time  interval  defined  by 

T < T < T 
1 2 


then  Eq.  (2.1)  gives 


= ar 2 

t -t  ex 

+ r 2 (b  - Si) 

l-ex 

2 i 

\ a / 

+ exf(T  ) 


(2.9) 

where 


u(t  ,r)  - u(t  ,r) 
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C(T  - T ) 

x = i 1— 

T 


At  late  times  when  the  displacement  becomes  constant, 
t can  be  much  larger  than  t so  that 

2 i 

ex-0 


a -*■  0 


b = u (°°,  r) 

The  steady-state  value  of  the  RDP  becomes 
n°°)  * r2  u (°°, r) 

This  result  also  follows  directly  from  Eq.  (2.1) 


(2.10) 


since 


?(•)  - 0 


(2.11) 


In  an  entirely  analogous  manner,  the  reduced  velocity 
potential  (RVP)  may  be  computed  from  the  velocity-time 
history  of  a particle.  If  u(t,r)  is  radial  particle 
velocity,  then 


u(t,r)  = - |-  llll  = iClI  + i£rl 
9r  r r2  ar 


(2.12) 


where  ’I'(t)  is  the  RVP. 

If  nr)  replaces  'F(t)  and  u(T,r)  replaces 
u(t,r)  in  Eq.  (2.9),  then  the  algorithm  for  calculating 
the  RVP  follows  directly. 

It  is  highly  desirable  to  transform  the  equivalent 
source  into  the  frequency  domain  in  order  to  identify  the 
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frequency  dependence  of  the  spectrum  required  for  yield 
scaling  (Eqs.  (2.5),  (2.7)  and  (2.8)).  Since  iW  , 0, 
then  T IS  bounded  at  zero  frequency  and  the  t-ourier  trans- 
orm  of  y(T)  may  be  done  numerically. 

Of  tb  A"  '"POrtant  relation  is  the  zero-frequency  limit 
f the  RVP  is  equal  to  the  steady  state  value  of  the  RDP 
Since 


V(f-O) 


00  00 
= / nr)  dr  = f 

J0 


«(T) 

dx 


dx 


Therefore,  since  ’{'(0)  = o, 


no)  = y(oo) 


(2.13) 


Reduced  velocity  potentials  have  been  computed  for  a 
variety  of  near  source  geologic  environments.  These  calcula- 

:™er\Perf0rn,ed  °n  3 Lagra"8ian  difference  code, 

ER,  that  simulates  a propagating  stress  wave  in  one 

space  dimension.  The  difference  equations  in  the  code  are 
identical  to  those  given  by  Cherry  and  Petersen  (1970).  The 
code  is  capable  of  carrying  the  propagating  stress  wave  into 
the  small  displacement,  elastic  regime  and  yet  flexible  enoug 
to  permit  appropriate  material  response  formulations  in  the 
large  displacement,  nonlinear  regime. 

Figure  2.1  shows  $ versus  frequency  for  a nuclear 
exp  osion  of  20  kT.  Each  source  function  in  the  figure  is 
i enti  led  by  an  integer  and  corresponds  to  the  20  kT  explo- 
sion detonated  in  a geologic  environment  with  a specific  set 
of  material  properties.  For  all  practical  purposes  the  spect, 
are  flat  (n=0)  within  the  teleseismic  frequency  band  (f  < J H; 
or  yields  up  to  at  least  200  kT.  This  implies  that  teleseis- 
mic body  wave  and  surface  wave  amplitudes  should  be  directly 
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proportional  to  the  yield  of  the  explosion.  Deviation  from 
this  linear  scaling  should  occur  within  the  yield  interval 
200  V T < W < 1,000  kT,  depending  on  the  source  function  used 
to  represent  the  rock  material  near  the  source. 

2-4  description  OF  THE  MATERIAL  MODFI. 

2.4.1  Introduction 

Our  goal  was  to  establish  both  the  shape  and  amplitude 
of  the  RVP  spectrum  as  a function  of  near  source  material 
properties.  The  properties  chosen  were 

f = mass  fraction  of  water  in  the  rock. 

<Pg  = the  initial  air-filled  volume  fraction. 

Eg  = the  elastic  pressure. 

Pc  = the  crush  pressure. 

kQ  *=  the  zero  pressure  bulk  modulus. 

P^  = the  overburden  pressure. 

P = the  shear  modulus. 

^m  ~ t^e  maximum  value  of  material  strength. 

Pm  = the  pressure  at  which  Ym  is  attained. 

em  = the  intetnal  energy  required  to  melt  the  rock- 
water  mixture. 

In  order  to  uncover  the  dependence  of  the  RVP  spectrum 
on  these  parameters  it  is  important  that  the  pressure  component 
P,  of  the  equation-of-state  and  the  material  strength,  Y,  be 
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uniquely  determined  when  these  parameters  are  specified. 
We  have  chosen  the  rock  constituent  to  be  NTS  tuff.  The 
following  sections,  2.4.1  and  2.4.2,  present  the  pressure 
and  material  strength  formulation  used  in  the  parameter 
study. 


2.4.2  The  Pressure  Component 


Given  the  pressure  component  of  the  equation  of  state 
for  both  water  and  grain  density  tuff,  then  the  pressure 
component  for  a fully  saturated  mix  is  obtained  by  specify- 
ing 


r. 


M 

f = „ — ~ - = mass  fraction  of  water  (2-14) 

r w 


For  water  we  used  the  EOS  formulation  of  Bjork  and 
Gittings  (1972)  while  for  tuff  we  used  the  formulation 
given  by  Riney,  £t  al . (1972)  with  a grain  density  of  2.4 
g/cc.  Figures  2.2  and  2.3  show  the  Hugoniot  and  release 
isentropes  for  these  two  constituents. 

At  a given  pressure  P a pressure  equilibrium  mix 
is  obtained  if 

v + f'v 
- _ r w 


e + f 'e 
r w 

e — rrr 


where 


V = 


r M 
f _ w 


rr  m 


(2.16) 


(2.17) 
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vr>  Vw’  "r  and  ew  are  sPecific  volumes  and  energies  of 
the  rock  and  water  components  at  the  pressure  F.  It  is  now 
possible  to  generate  a taole  that  gives  P as  a function  of 
v and  e.  This  table  (Riney,  et  al^  (1972))  suitably  chooses 
values  of  v and  e so  that  no  time  consuming  search  is 
needed  when  performing  an  interpolation  at  a given  thermo- 
dynamic state.  Figure  2.4  shows  the  result  of  mixing  the 
rock  and  water  constituents  shown  in  Figs.  2.2  and  2.3  using 
a 0.17  value  for  f. 

Air-filled  porosity  is  included  in  the  model  by  de- 
fining a parameter  a,  given  by 


a 


V 

r 


+ V + V 
w 

V 

r w 


(2.18) 


where  v is  the  specific  volume  of  the  partially  saturated 
mix.  Then 


p = | P(v/a,  e) 


(2.19) 


where  P is  found  from  the  table  using  v,  a,  e.  Equation 
(2.19)  is  equivalent  to  the  porosity  model  proposed  by 
Herrmann  (1969)  except  F is  weighted  with  1/a.  Both  v 
and  e are  zone  variables  and  are  obtained  from  SKIPPER  at 
any  given  time  cycle. 


The  parameter  a is  assumed  to  vary  in  the  following 
way  during  loading  (pore  collapse)  (Cherry,  et  al.,  1972) 


a = 1 P > P 

c 


(2.20a) 


a = 1 + (a 
v e 


1)  1 


P-P  \2 

H H 

c e/ 


P < P < P 
e - - c 


(2.20b) 
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/ n p/p  \ 

“ = “o  + B l1  ' e ) 0 1 p 1 pe  (2.20c) 

a = PaP=0  + a0  P < 0 (<*'  = da/dp)  (2 . 20d) 

where  Pg  and  Pc  are  the  elastic  and  crush  pressure.  B 

and  n are  determined  by  specifying  the  zero  pressure  bulk 
modulus  (k^)  of  the  porous  mix  and  requiring  that  a'  be 
continuous  at  P - Pg.  The  initial  air-filled  volume  frac- 
tion (4>  ) and  a are  related  bv 
o o ' 

a - 1 V 

+ = —2 - p 

o a V + V + V C 2 . 21) 

c w r p 

while  the  volume  fractions  of  rock  and  water  (4>  and  <f>  ) 
are  given  by  r w 

_ 1 - ♦ 

*r  " 1 + A°  C 2 . 2 2) 

^w  = (l  - A°  ) A C2.23) 


where 


p . Pr  = 2.4  g/cc 

A = -I  * 

p«  1_  ”w  - 1-0  g/cc 

We  have  found  (Cherry,  et  aL  (1972))  that  a pore  re- 
covery formulation  needed  to  be  included  in  the  model  in  order 
to  insure  a smooth  transition  between  loading  states  near  P 
If  a*  is  the  minimum  a seen  by  a zone  then 

a = ta*  + (1  - t)  a (P)  n 
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where  a^P)  is  given  by  Eqs.  (2.21)  and 

t = 0,  a*  > ag  = -y  - e , 1 < a*  < a (2.25) 

e c 

C 

This  assures  full  recovery  when  a*  = a and  zero  recovery 
when  a*  = 1 . 

Figure  2.5  shows  the  loading  and  release  P-v  curves 
that  result  from  adding  5 percent  air-filled  porosity  to  the 
mix  of  Fig.  2.4.  We  assume  Pg  = 0.15  kbar  and  P = 1.25 
kbar.  Also  shown  for  comparison  in  the  figure  are^he  hydro- 
static unloading  data  for  unit  3 tuff  at  the  Diamond  Dust 
site  (Stephens,  et  al^  (1970)). 

In  summary,  the  pressure  component  of  the  equation  of 
state  is  uniquely  determined  by  assuming  a pressure  equili- 
brium mix  and  specifying  f (Eq.  (2.14)),  <p  (Eq.  (2.21)), 

^e*  V k0  (2.20b)  and  (2.20c))  and  the  grain  density 

of  the  rock  component . 


2.4.3  Material  Strength 

A strength  relation  of  the  form 


SijSij  < f V-  (PT.e5  (2.26 

where  S..  is  the  stress  deviator,  PT  = P + P , p equals 
the  overburden  pressure,  and 


S; 


Y(P,e)  = p 

pt1 

2 - r 

1 - ®-‘ 
e 

m 

L m 

P < P 

m 


P > P 
— m 


= 0 


e > e 
— m 


(2.27) 
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C por  spherical  symmetry  Eq.  (4.26)  becomes 


F 


i 


I ^ 


E 


I 0 


5 

f 


I 


f 


i- 


|s, I £ § Y (PT,e) 


(2.28) 


For  all  the  calculations  discussed  here,  a generalized 
associated  flow  rule  was  used  for  the  material  in  which  the 
volumetric  plastic  strain  rate  was  assumed  to  be  zero.  This 
assumption  implies  (Cameron  and  Scorgie  (1968))  that  if 
|Sjl>  calculated  from  Eq.  (2.20),  exceeds  | Y,  then 

Si  = (f  Y)  sign  (2.29) 

The  material  strength  dependence  on  P , Y and  P 
is  summarized  in  Fig.  2.6.  The  shock  loading  path  is  also 
shown  in  the  figure  to  have  a slope  of  y/k,  where  k is 
the  effective  bulk  modulus  of  the  shock  wave.  This  loading 
path  is  only  valid  for  regions  where  the  radial  strain  is 
much  greater  than  the  tangential  strain,  i.e.,  at  the  shock 
front. 


2- 5 results  of  the  parameter  study 
2.5.1  Introduction 


The  results  presented  in  Fig.  2.1  show  the  calculated 
RVP  spectra  using  the  material  model  presented  in  Section 
2.4.  As  noted  earlier,  the  spectra  are  essentially  flat  in 
the  teleseismic  frequency  band  for  device  yields  at  least 
as  large  as  200  kT.  This  result  along  with  Eq.  (2.4),  (2.4), 
and  (2.13)  imply  that  surface  wave  magnitude,  m^,  has  the 
proportionality 


Z 


► 


mb  ^ log 

while  surface  wave  magnitude,  Mg,  obeys 
Ms  * log  [y  ^(oo)] 


11=2 

a 


(2.30) 


(2.31) 
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Y * Maximum  allowable  stress  difference 

PT  = hydrodynamic  component  of  the  stress  tensor 
including  overburden  pressure 

Pq  = Overburden  pressure 


Figure  2.6. 


Assumed  relationship 
strength  (Y)  and  the 
of  stress  C?T) • 


between  the  material 
hydrodynamic  component 
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If  4'iC“)  and  Yk(»)  are  two  values  of  Y(«)  cor- 
responding to  materials  i and  k then  the  change  in  tele- 
seismic  magnitude  is  given  by 


Am 


m 


- ni  = log 


V08) 


where 


ak 

a - - for  mb 


(2.32) 


If  material  properties  remain  constant  and  the  only  change 
m device  yield  then  according  to  Eqs.  (2.7)  and  (2.8) 

i V W. 

Am  = m - in  = log  -1 

Wk 

Equation  (2.32)  gives  the  "scaling"  law  for  teleseis- 
mic  magnitudes  in  terms  of  ¥(«.).  This  equation  is  always 
valid  as  long  as  the  spectrum  of  the  equivalent  source  is 
flat  within  the  teleseismic  frequency  band  and  path  effects 
associated  with  earth  structure  are  invariant  between  events. 

As  an  example,  assume  that  the  values  of  y(«)  for 
a yield  of  0.02  kT  in  tuff  and  granite  are  8 m3  and  12.5  m3 
with  corresponding  values  of  a being  2.6  km/sec  and  5.6 
km/sec,  and  y of  40  kbar  and  275  kbar  respectively,  then 

”b  • mb  ■ l0S  276  • l0«  Htf  = °-14 


I 
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while 


MS  ' Ms  “ lQg  C40)(8)  - log  (275) (12 . 5)  = -1.03 

While  both  and  Mg  are  measures  of  T(oo)  they  also 
are  severely  affected  by  a and  y.  The  above  example 
serves  to  illustrate  the  importance  of  these  near  source 
elastic  parameters  on  seismic  coupling. 

2.5.2  The  Sensitivity  of  ¥(*)  to  Air-Filled  Porosity, 
Material  Strength  and  Overburden  Pressure 

Table  2.1  lists  the  material  properties  (d>  , P , p 

0 G J c 9 

Ym’  V Po}  along  with  the  calculated  values  of  V(°°)  for  an 
explosion  yield  of  0.02  kT.  Also  included  in  the  table  is 
Am,  computed  from  Eq.  (2.32)  where  calculation  5 has  been 
chosen  for  the  reference  magnitude.  The  value  of  "a"  in 
Eq.  (2.32)  is  unity  (ot  = 2.4  km/sec,  y = 40  kbar  for  all 
calculations)  and  f = 0.17. 

The  following  expansion  summarizes  the  results  of  the 
parameter  study 

Am  = (-20.2  + 355  <p  ) A<|>  - 0.38  Y (kbar) 

oo  m ' 

- 0.34  P'1  AP  + 2.0  A?  (kbar)  + 0.057  AP  (kbar) 
oo  e c ' 

+ 0.41  APm  (kbar)  (2.33) 

Increasing  <p^,  and  P^  decreases  teleseismic  magnitude 
while  increasing  Pg,  P^  and  increases  teleseismic  magni- 
tude . 

Figure  2.7  shows  Am  versus  <p  for  calculations  1, 

6 and  4.  Am  is  highly  sensitive  to  change  in  4>  . Over 
the  internal  0 < $ < 0,02,  Am  decreases  by  0.313  units 
and  the  seismic  coupling  is  reduced  by  a factor  of  two. 


SUMMARY  OF  SOURCE 
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Figure  2.8  shows  the  influence  that  the  maximum  al- 
lowable material  strength,  Y , has  on  magnitude.  Increasing 
^m  ^ a factor  of  two  causes  magnitude  to  decrease  by  0.102 
magnitude  units.  Therefore,  the  device  yield  could  be  in- 
creased by  a factor  of  1.55  in  the  strong  material. 

These  calculations  show  that  a porous,  high  strength 
material  optimizes  teleseismic  decoupling.  Increasing  both 
the  porosity  and  material  strength  may  obscure  the  device 
yield  by  at  least  a factor  of  three  when  teleseismic  magni- 
tube  is  used  to  infer  yield.  We  should  note  that  both  high 
air-filled  porosity  and  high  material  strength  imply  that 
the  material  characterizing  the  near  source  environment  is 
above  the  water  table. 

For  a given  value  of  Y , the  shape  of  the  material 
strength  surface  is  controlled  by  P^,  the  pressure  at  which 
Ym  is  attained*  From  calculations  9 and  7 we  see  that  in- 
creasing Pm  by  a factor  of  2.08  causes  magnitude  to  in- 
crease by  0.26  units.  The  coupling  efficiency  increases  as 
Pm  ls  increased,  which  is  opposite  to  the  effect  of  in- 
creasing or  Ym.  This  simply  indicates  that  the  shape 

of  the  failure  surface  plays  an  important  role  in  determining 
coupling  efficiency.  The  material  has  less  strength  at  low 
pressures  when  Pm  is  increased  and  we  should  expect  a cor- 
responding increase  in  magnitude. 

The  effect  of  overburden  pressure,  P , on  magnitude 
is  shown  in  Fig.  2.9.  Coupling  efficiency  decreases  with 
increasing  overburden  pressure,  when  all  other  material 
properties  remain  constant.  An  increase  in  the  depth  of 
burial  (DOB)  by  a factor  of  six  causes  magnitude  to  decrease 
by  0.3  units,  which  in  turn  permits  a yield  increase  by  a 
factor  of  two.  The  decrease  in  magnitude  is  not  uniform  with 
depth.  An  increase  in  DOB  at  shallow  depths  causes  magnitude 
to  decrease  more  then  that  caused  by  an  equivalent  increase 


allowable  material  strength  on  magnitude. 


in  DOB  a*  greater  depths.  This  effect  results  from  the  rapid 
increase  in  material  strength  at  low  pressures. 

We  should  note  that  it  is  unlikely  that  material  pro- 
perties would  remain  invariant  as  overburden  pressure  is  in- 
creased from  100  bars  to  600  bars.  An  increased  DOB 
usually  corresponds  to  an  increase  in  water  content,  result- 
ing in  decreased  air-filled  porosity  and  material  strength. 
These  changes  would  more  than  compensate  for  the  overburden 
pressure  decoupling  effect. 

Calculations  9 and  10  show  that  magnitude  is  not  very 
sensitive  to  large  changes  in  P , the  pressure  at  which  all 
ambient  porosity  is  irreversibly  removed.  An  increase  in 
Pc  by  a factor  of  2.5  produces  an  increase  in  magnitude  of 

0.03  units. 

Finally,  calculations  4 and  3 show  the  sensitivity  of 
magnitude  to  P , the  pressure  15  lit  at  which  irreversible 

c 

pore  collapse  begins.  An  increase  in  Pg  produces  an  in- 
crease in  coupling  efficiency.  A preshot  estimate  of  Pg 
that  is  incorrect  by  a factor  of  two  will  result  in  an  error 
in  predicted  magnitude  by  0.2  units. 

2.6  SUMMARY  AND  COMMENTS 

1.  All  RVP  spectra  calculated  to  date  have  been  flat 
in  the  teleseismic  frequency  bend  for  device  yield 
at  least  as  large  as  200  kT.  For  this  yield  inter 
val,  both  body  wave  and  surface  wave  amplitudes 
should  scale  linearly  with  device  yield. 

2.  The  near  source  elastic  constants  which  influence 
teleseismic  amplitudes  are  the  P wave  velocity,  a, 
and  the  shear  modulus,  y.  Body  wave  amplitude  is 
inversely  proportional  to  a.  Surface  wave 
amplitude  is  directly  proportional  to  y. 
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Increases  in  air-filled  porosity  ($  ),  maximum 
material  strength  (Ym)  and  overburden  pressure 
(P  ) causes  a reduction  in  seismic  coupling. 
Equation  (2.33)  shows  that  increasing  <J>  from 
0 to  0.01  causes  magnitude  to  decrease  by  0.2 
units.  Increasing  by  0.5  ;bar  also  causes 

magnitude  to  decrease  by  0.2  units.  Magnitude 
dependence  on  overburden  pressure  is  given  by 
(-0.34  AP^/P^).  At  shallow  burial  depths  magni- 
tude is  very  sensitive  to  small  changes  in  over- 
burden pressure.  For  example,  if  P = 0.075  kbar 
and  AP^  = 0.044  kbar  then  magnitude  decreases  by 
0.2  units.  This  change  in  overburden  pressure 
corresponds  to  a change  in  DOB  of  approximately 
650  feet. 

Increases  in  the  elastic  pressure  (P  ),  the 
crush  pressure  (Pc)  and  the  pressure  at  which  the 
material  reaches  its  maximum  strength  (P  ) cause 
an  increase  in  seismic  coupling.  Equation  (2.33) 
shows  that  increasing  Pg  by  0.1  kbar  causes  a 
magnitude  increase  of  0.2  units.  Increasing  P 
by  0.5  kbar  also  causes  a magnitude  increase  of 
0.2  units.  Magnitude  is  not  very  sensitive  to 
changes  in  P . 

The  results  summarized  by  Eq . (2.33)  are  dependent 
on  the  material  model  assumed  for  the  rock-water 
mixture.  Neither  strain  rate  effects  on  brittle 
fracture  of  the  material  as  the  failure  surface 
is  exceeded  have  been  included.  The  latter  is 
certain  to  be  important  at  shallow  depths  of 
burial.  The  former  may  play  a major  role  in 
broadening  the  pulse  shape  and  reducing  the  yield 
range  over  which  the  spectrum  is  flat. 
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III.  A THEORETICAL  INVESTIGATION  OF  THE  TELESEISM I C 
P WAVE  FROM  UNDERGROUND  EXPLOSIONS 

3.1  INTRODUCTION 

Using  our  deterministic  computational  techniques,  we  have 
been  attempting  to  identify  and  understand  the  important  features 
controlling  the  teleseismic  signature  of  underground  explosions. 
For  short  period  P waves,  some  interesting  results  obtained  dur- 
ing the  last  year  are  summarized  below. 

In  Section  3.2  the  influence  of  the  source  spectrum  on 
the  shape  (frequency  content)  of  the  teleseismic  P wave  signa- 
ture of  explosions  is  discussed.  We  also  discuss  the  dependence 
of  amplitude-yield  scaling  on  various  source,  propagation  path 
and  measurement  techniques. 

Some  interesting,  and  perhaps  paradoxical  consequences 
of  the  standard  formulation  of  the  anelastic  attenuation  of 
teleseismic  waves  are  discussed  in  Section  3.3. 

In  Appendix  E is  given  a theory  whereby  an  arbitrary 
explosion  or  earthquake  source  may  be  embedded  in  a plane  layered 
model  of  the  earth's  crust.  This  theory  is  applied  in  Section 
3.4  to  investigate  the  influence  of  near  source  crustal  struc- 
ture on  the  teleseismic  P wave.  Finally,  in  Section  3.5  we 
study  the  influence  of  tectonic  release  superposed  on  the 
spherically  symmetric  explosion  source. 

3.2  THE  DEPENDENCE  OF  THE  TELESEISMIC  P WAVE  TRAIN  ON  THE 

SOURCE  SPECTRUM  FOR  EXPLOSIONS  IN  NTS  TUFF 

We  know  that  the  teleseismic  P wave  signatures  of  under- 
ground explosions  are  affected  by  the  characteristics  of  the 
source  and  its  immediate  vicinity,  the  travel  path  through  the 
upper  mantle,  the  geology  of  the  receiver  vicinity,  etc. 

When  possible,  one  tries  to  isolate  the  various  important  para- 
meters for  separate  study.  In  this  section  the  intention  is 
to  study  the  effect  of  the  source  and  its  immediate  vicinity. 
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A set  of  observed  seismograms  for  five  explosions 
located  within  a few  kilometers  of  each  other  in  the  same  re- 
gion at  NTS  is  shown  in  Fig.  3.1.  The  yields  of  these  five 
explosions  vary  by  at  most  a factor  of  two  while  the  burial 
depths  are  the  same  to  within  10  percent.  All  five  recordings 
are  from  an  array  of  short-period  sensors  located  more  than 
30  degrees  from  NTS.  While  the  recordings  are  similar,  they 
have  distinct  differences.  Due  to  the  similarities  of  the 
events,  these  differences  can  only  be  accounted  for  by  the 
source  itself  and  its  immediate  vicinity  --  to  a depth  of  one 
or  two  kilometers.  We  first  examine  the  possible  influence  of 
the  ’.ource  spectrum  on  the  observed  signals. 

The  results  of  a series  of  one-dimensional  explosion 
source  calculations  for  material  parameters  typical  of  NTS  tuff 
have  been  described  in  previous  S3  reports  (e.g.,  Cherry,  et  al. 
(1974b),  Bache , ejt  al . , (1974)).  As  was  pointed  out  in  those 
reports,  the  far  field  component  of  the  source  generated  P wave 
is 

77  ( d _ (u)  _-iwR/a 

US(R’u)  ~ “fcr  e (3.2) 

where  a is  P wave  velocity  and  R is  the  radial  coordinate. 
The  source  spectrum  is  then  given  by  ¥(0)),  the  reduced 
velocity  potential  transform.  This  quantity  is  plotted  for 
a range  of  material  properties  typical  of  NTS  tuff  in  Fig. 

3.2.  The  calculations  portrayed  in  the  figure  are  those 
of  the  parameter  study  described  in  the  above  named  references. 

Most  apparent  in  Fig.  3.2a  is  the  essential  similarity 
of  the  amplitude  spectrum  for  a wide  range  of  coupling  ef- 
ficiencies. That  is,  while  the  zero  frequency  limit,  ¥ , 
varies  by  as  much  as  a factor  of  5.5,  the  basic  shape  of 
the  spectra  are  relatively  unchanged. 

In  Fig.  3.2b  a detailed  look  at  one  of  these  spectra, 

#10  if  3.2a,  is  given  on  a log-log  plot,  along  with  an 
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indication  of  the  rate  of  falloff  beyond  the  corner.  All 
the  source  calculations  were  carried  out  at  the  same  yield. 
Cube  root  scaling  is  applied  to  obtain  source  spectra  at 
differing  yields. 

The  source  #10  of  Fig.  3.2  is  scaled  to  a wide  range 
of  yields  from  10  to  4,000  kT  in  Fig.  3.3a.  It  is  seen 
that,  for  the  lower  yields,  the  source  spectrum  is  flat  over 
the  entire  band  of  frequencies  (0.5  - 2.0  Hz)  of  interest  for 
short  period  P waves.  For  high  yields  the  corner  does  ap- 
pear in  the  teleseismic  frequency  band.  However,  attenuation 
by  Q tends  to  obscure  this  effect. 

In  Fig.  3.3b  the  source  spectra  of  3.3a  are  multiplied 
by  exp  [ -irf ] . That  is,  T/Q  is  taken  to  be  unity.  This 
value  of  T/Q  has  been  found  by  our  studies  and  by  others 
(e.g.,  Helmberger  (1973),  Carpenter  (1966))  to  be  reason- 
able for  most  teleseismic  P wave  phases. 

The  direct  P waves  represented  by  the  spectra  of 
Fig.  3.3  are,  of  course,  seen  by  an  instrument  which  further 
masks  details  of  the  source  spectra.  In  Fig.  3.4a  the 
spectra  of  3.3b  are  multiplied  by  the  instrument  transfer 
function  for  an  LRSM  short  period  seismograph.  As  seen 
through  the  filter  of  the  instrument  and  the  anelastic 
earth,  the  source  amplitude  spectra  are  almost  identical  in 
shape  over  this  broad  yield  range.  To  provide  an  upper 
bound  on  the  kinds  of  source  amplitude  spectral  differences 
expected  from  yield  scaling  alone,  the  same  information  is 
plotted  in  Fig.  3.4b,  this  time  for  T/Q  = 0.5.  Values  of 
Q this  high  are  not  expected  to  be  encountered  for  teleseis- 
mic P waves  from  NTS  explosions. 

It  should  be  emphasized  that  the  spectra  shown  thus 
far  are  for  the  direct  P wave  only.  That  is,  modification 
by  reflection  from  the  free  surface  and  other  near  source 
velocity  discontinuities  is  not  included.  The  effect  of  the 
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most  elementary  modification  by  reflected  phases  is  shown 
in  Fig.  3.5.  First,  the  source  spectra  for  the  six  yields 
of  Fig.  3.3a  at  scaled  depths  of  burial  arc  show  including 
the  direct  and  free  surface  reflected  (pP)  phases.  The  re- 
flection coefficient  was  taken  to  be  0.9  which  is  typical 
for  takeoff  angles  for  telcseismic  P waves.  The  burial 
depths  were  cube  root  scaled  to  an  assumed  depth  of  300  m 
for  the  10  kT  explosion.  The  modification  of  these  source 
spectra  by  the  LRSM  instrument  and  Q (T/Q  = 1.0)  is  indicated 
in  Fig.  3.5b. 

These  amplitude  spectral  density  plots  clearly  indi- 
cate how  source  spectrum  details,  such  as  corner  frequency 
or  rate  of  amplitude  decay  beyond  the  corner  for  teleseismic 
P wave  spectra,  are  almost  totally  obscured  by  the  filter 
of  the  earth  and  instrument  (in  principle,  both  can  be  re- 
moved) and,  more  importantly,  by  the  presence  of  later 
arriving  phases  such  as  pP. 

Source  spectra  having  a pronounced  hump  near  the 
corner  (#4  of  Fig.  3.2  exhibits  this  behavior  to  a small 
degree)  have  been  suggested  by  some  investigators  (e.g., 

Werth  and  Herbst  (1963)).  Presence  of  such  an  overshoot 
of  moderate  size  would  have  negligible  effect  on  the  results 
presented  above.  Aki,  £t  al . (1974)  have  inferred  the  pre- 
sence of  an  overshoot  of  as  much  as  five  times  the  low  fre- 
quency value  and  tentatively  attributed  this  to  spallation, 
an  effect  not  included  in  our  calculations.  A huge  over- 
shoot of  this  kind,  should  it  cube-root  scale  with  yield, 
would  have  a noticeable  effect  on  the  teleseismic  spectrum. 

Examination  of  amplitude  spectral  density  plots  are 
useful,  but  can  be  misleading  since  the  phase  component  of  • 
the  spectrum  is  ignored.  Therefore,  we  move  on  to  an  examina- 
tion of  results  in  the  time  domain;  that  is,  synthetic 
seismograms . 
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' 1 lustrative  synthetic  seismograms  are  given  in  Figs. 
3.6  through  3.9.  For  each  seismogram  the  source,  #10  of 
Fig.  3.2,  was  embedded  at  the  indicated  depth  of  burial  in 
a crustal  structure  typical  of  NTS.  The  theory  of  Appendix 
B was  used  with  the  crustal  structure  of  Table  3.1. 


TABLE  3.1 


NTS  CRUSTAL  STRUCTURE  FOR  THE  SYNTHETIC  SEISMOGRAMS 
OF  FIGURES  3.6  THROUGH  3.9 


Depth 
2.5  km 
18.0 
31.0 


P Velocity 
2.4  km/sec 
6.0 
6.7 
7.9 


S Velocity 
1.47  km/sec 

3.5 
3.85 

4.6 


Density 
1.9  g/cm3 
2.8 
3.0 
3.3 


The  seismograms  are  for  the  first  mantle  arrival  at 
a distance  of  4,000  km.  A geometric  spreading  factor  of 
0.63  x 10- 4 for  the  travel  path  between  the  base  of  the 
crust  at  the  source  and  the  earth's  surface  it  the  receiver 
was  computed  by  ray  theory  and  applied.  The  Q operator 
including  dispersion  as  formulated  by  Strick  (1970)  and  dis- 
cussed in  Section  3.3  was  applied. 

In  Figs.  3.6  and  3.7  the  effect  of  yield  scaling  of 
the  source  spectrum  is  shown  by  seismograms  at  a wide  yield 
range  but  at  a single  depth  of  burial.  For  these  seismograms 
we  have  taken  T/Q  = 1.1.  The  seismograms  of  Fig.  3.7  are 
repeated  m Fig.  3.8  for  the  case  of  a Q which  is  twice  as 
high  (T/Q  = 0.55).  Finally,  in  Fig.  3.9  the  seismograms  at 
scaled  depth  of  burial  are  shown.  Except  for  the  somewhat 

lower  Q,  the  latter  set  correspond  directly  to  the  spectra 
of  Fig.  3.5. 
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The  synthetic  seismograms  shown  in  Figs.  3.6  through 
39  were  used  to  generate  amplitude-yield  scaling  relation- 
ships. These  seismograms,  which  correspond  to  four  different 
p ysical  situations  previously  described,  allow  us  to  examine 
the  dependence  of  amplitude-yield  scaling  on  factors  such 
as  depth-of-burial,  anelastic  attenuation  (Q)  and  the  details 

of  the  manner  in  which  the  actual  amplitude/period  measure- 
ments  are  made. 

Amplitude  and  half-period  measurements  were  made  for 
three  different  wavelets  Ca„,  b and  c)  for  each  of  the  P 
wave  trains  in  Figs.  3.6  through  3.9.  All  measurements  were 
taken  from  computer  printouts  of  the  seismograms,  rather  than 
e p Ots  in  Figs.  3.6  to  3.9,  for  increased  precision.  The 
e lnitions  of  aQ,  b and  c (Fig.  3.6a)  are  the  same  as  those 
given  by  Springer  and  Hannon  (1973).  In  addition  to  these 

measurements  the  rise-times  of  the  a wavelet  were  also 
determined. 

tnde  e"  Fl8!’  3‘10’  3-11  a"d  3,12  we  have  Plotted  the  ampli- 
of  a0,  b and  c versus  yield,  respectively.  For  any 

given  figure  there  are  three  different  types  of  symbols  plot- 

e ■>  A » and  * )•  The  convention  adopted  for  these 

symbols  is  the  following:  . represents  amplitude  measurements 

uncorrected  for  instrument  response;  a represents  amplitudes 

lvided  by  the  instrument  response  factor  appropriate  to 

twice  the  half-period  of  the  particular  wavelet  under  consider.- 

tion,  the  . represents  amplitudes  divided  by  the  product  of 

he  instrument  response  factor  and  twice  the  half-period 

of  the  wavelet,  had.  figure  in  the  sequence  3.10,  3.11,  3.12 

four  C 5 corresP°uding  to  the 

3s:rions  of  de,,th-6f-‘°uriai  ^ —ib.d 

Considering  the  a0  wavelet  of  Figs.  3.10a-d,  there  are 
several  points  to  be  made.  First,  for  the  seismograms  of 

lg5'  3'7  ' 3'8  corresponding  to  3.10b-c,  the  maximum  duration 


52 


(°»)  apnjixdiuv 


1 

• • • t • * 

. 

nrit: 

[Jj  . 

1L 

I 

l 

in  , 
;.!i 

i 

i 

. ! J • 

’t 

• . » 1 - 

. 

i 

■ r 

jliil 

Us. 

/ 


f 


J 


scaled  DOB 


mm 


immmppp 


('•!)  J|>ni  l 


rttfh 


(q)  apiuiiduiy 


<L»  PQ 


(q i ;diuv 


R- 2646 


(0.775  s*c)  of  this  wavelet  is  less  than  the  pP-P  delay 
time  (0.83  sec).  Thus,  these  two  cases  demonstrate  the  de- 
pendence of  amplitude-yield  scaling  on  details  of  the  mea- 
surements alone.  The  situation  in  Fig.  3.6  or  3.10a  is  just 
the  opposite;  that  is,  the  pP-P  delay  time  (0.42  sec)  is 
less  than  the  minimum  duration  (0.525  sec)  of  the  aQ  wavelet 
which  implies  that  aQ  is  a mixed  (pP+P)  phase.  In  the  case 
of  Fig.  3.9  or  3.10d  (scaled  depth-of-burial),  aQ  is  a mixed 
phase  for  yields  of  10  and  100  kT  but  not  for  the  higher 
yields . 

The  amplitudes  (■),  uncorrected  for  instrument  res- 
ponse and  dominant  period,  of  the  aQ  wavelet  in  Figs.  3.10a-c 
show  a marked  departure  from  a linear  dependence  on  yield  at 
yields  > 100  kT.  This  behavior,  particularly  for  the  situa- 
tions of  Figs.  3.10b  and  c where  depth-of-burial  is  not  in- 
fluencing the  results,  is  ascribed  to  the  explosion  corner 
frequencies  shifting  to  lower  frequencies  with  increasing 
yield  and  sliding  into  the  passband  of  the  instrument  (LRSM 
in  this  case).  This  explanation  has  been  invoked  by  others 
(Haskell,  (1967) , Evernden  and  Filson,  (1971))  for  similar  be- 
havior of  later  portions  of  the  P-wave  trains  from  explosions. 
Simple  correction  of  these  amplitudes  for  instrument  response 
at  the  dominant  period  (equal  to  twice  the  duration  of  a ) 
results  in  approximately  linear  amplitude-yield  relationships 
(note  the  trends  of  the  a's  in  Figs.  3.10a-c).  Further 
manipulation  of  the  amplitudes,  namely  division  by  the  domi- 
nant period,  reintroduces  a flattening  in  the  data  at  higher 
yields.  This  last  result  is  quite  significant  since  amplitude/ 
period  is  the  quantity  employed  in  conventional  mb  estimates 
and  the  data  in  Figs.  3.10a-c  indicate  that  even  in  the  absence 
of  such  effects  as  depth-of-burial,  mb  will  not  exhibit  a 
linear  dependence  on  explosion  yield.  The  quantity  that 
should  be  compared  with  explosion  yield  is  amplitude  corrected 
for  instrument  response;  that  is,  apparent  ground  motion. 
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Dio  scaling  relationships  shown  Cm  :i  (Fig.  | o ) 
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tivfly,  are  complicated  hy  depth-of  - bur  in  1 effects,  particu- 
larly in  Figs.  3.10d,  3. lid  and  3.12d  corresponding  to  the 
situation  where  the  explosions  are  detonated  at  scaled 
depths.  The  c wavelet  in  Figs.  3.10a  and  b corresponds  to 
the  maximum  amplitude  in  the  P-wave  trains  considered  here 
(see  Figs.  3.6  and  3.7)  and  would  be  used  for  conventional 
mb  estimates.  The  data  in  Figs.  3.12a  and  b suggest  that  the 
amplitude  of  ground  motion  rather  than  amplitude/period  once 
again  results  in  a more  nearly  linear  dependence  on  yield. 

The  synthetic  seismograms  shown  in  Figs.  3.8  (T/Q  = 
0.55,  DOB  = 1 km)  and  Fig.  3.9  (scaled  depth-of -bur ial) 
demonstrate  one  of  the  major  problems  in  magnitude  esti- 
mates for  explosions.  The  maximum  amplitudes  of  the  first 
few  cycles  in  these  figures  are  sometimes  associated  with  the 
b wavelet  and  other  times  with  the  c wavelet.  As  expected, 
this  variation  contaminates  the  amplitude-yield  relation- 
ships as  can  be  clearly  seen  in  Fig.  3. lid  and  Figs.  3.12c 
and  d.  While  the  amplitude  of  the  aQ  wavelet  is  the.  pre- 
ferred measurement  for  yield  estimates,  signal -to-noise 
problems  might  be  too  severe  in  the  case  of  low  yield  explo- 
sions recorded  at  teleseismic  distances.  This  being  the 
case,  then  the  data  shown  here  support  a preference  for  the 
b wavelet  rather  than  the  c or  maximum  cycle  in  the  P-wave 
train. 

As  previously  remarked,  the  duration  of  the  aQ  wave- 
let is  relatively  (with  few  exceptions)  uncontaminated  by 
such  effects  as  the  free  surface  reflection  (pP)  and  receiver 
crustal  structure.  As  a result, this  wavelet  was  chosen  for 
a special  investigation  of  the  behavior  of  dominant  period 
with  explosion  yield  as  modified  by  attenuation  (Q)  and 
instrument  response. 
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T*'o  duration  of  nn  was  determined  for  a series  of 
explosions  detonated  at  Pahute  Mesa.  The  se  i smog  rains  used 
were  recorded  at  one  of  the  short-period  stations  of  the 
Palmer  network  in  Alaska.  Because  of  classification  con- 
siderations a fit  (solid  line)  to  the  observed  data,  rather 
than  the  data  set  itself,  is  given  in  Fig.  3.13.  The  verti- 
cal bars  with  arrows  enclosed  give  a measure  of  the  precision 
of  these  data.  The  one-sided  arrow  with  AtQ  above  is  an 
estimate  of  the  maximum  amount  by  which  the  duration  of  the 
a0  wavelet  may  be  underestimated  for  the  lower  yield  explo- 
sions studied.  This  underestimate  arises  because  of  the 
tendency  to  pick  a late  arrival-time  for  the  onset  of  the 
a0  wavelet  in  a background  of  seismic  ncise.  While  the  error 
is  on  the  order  of  0.1  sec  for  explosions  of  100  kT  yield,  it 
is  considerably  less  (<  0.05  sec)  for  higher  yield  (>  1 MT) 
explosions.  The  dashed  line  indicates  a "fit"  to  the  data 
with  AtQ  added  to  estimates  of  the  duration  of  a0,  thus 
illustrating  the  bound  on  the  possible  range  of  the  dura- 
tion underestimates.  The  remaining  points  on  Fig.  3.13 
(see  legend)  are  from  theoretical  seismograms  and  corres- 
pond to  the  four  different  combinations  of  depth  of  burial 
and  attenuation  (T/Q)  described  previously. 

While  not  conclusive,  there  are  a few  points  that  we 
can  make  about  the  results  of  Fig.  3.13.  First,  the  pre- 
dicted behavior  for  the  case  of  T/Q  = 0.55,  DOB  = 1 km  (0) 
is  definitely  at  variance  with  the  observations.  However, 
for  the  same  DOB  (1  km)  but  T/Q  = 1.1  the  predicted  points 
(o)  actually  follow  the  observed  trend  within  the  uncer- 
tainty of  the  observations.  In  the  case  of  scaled  depth 
of  burial  (□)  the  point  at  100  kT  plots  outside  the  range 
of  observations.  While  the  actual  Pahute  Mesa  explosions 
were  approximately  buried  at  scaled  depths  for  containment 
purposes,  it  has  been  reported  (Frasier  (1972);  Springer 
(1974))  that  the  pP-P  delay  times  for  Pahute  Mesa  events 
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Yield  versus  duration  of  ag  wavelet  for  pre- 
dicted and  observed  seismograms.  Measurement 
for  Q and  □ not  available  at  1000  kT. 
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are  considerably  greater  than  those  used  in  the  computations. 
If  this  is  correct,  then  the  open  circle  turns  out  to  be 
more  appropriate,  since  for  that  case  the  pP-P  delay  is 
greater  than  the  duration  of  a^,  whereas  in  the  case  of  the 
100  kT  explosion  at  a scaled  depth  of  burial  the  a^  wavelet 
was  contaminated  by  pP. 

An  investigation  such  as  the  one  described  in  this 
section  is  quite  instructive  for  pointing  out  optimum  pro- 
cedures for  yield  estimation  of  underground  explosions. 

It  allows  one  to  isolate  those  factors  that  most  strongly 
affect  the  frequency  content  and  amplitude  of  observed 
explosion  signatures  and  thereby  correct  amplitude-yield 
scaling  relationships  accordingly. 

3.3  THE  DISPERSION  ASSOCIATED  WITH  A CONSTANT  Q,  LINFAR 

ABSORPTION  OPERATOR 

It  is  universally  accepted  that  seismic  waves  are 
attenuated  anelastically  during  their  passage  through  the 
earth.  The  mechanism  of  this  attenuation,  and  even  its  form 
over  the  measurable  frequency  band  is,  however,  far  from 
being  fully  understood.  For  those  interested  in  the  con- 
struction of  synthetic  seismograms,  the  common  practice  is 
to  assume  that  the  attenuation  is  satisfactorily  represented 
by  a frequency  independent  Q and  apply  the  operator 
exp  [ - uiT / 2Q j where  T is  the  source -receiver  travel  path  and 
Q is  the  effective  attenuation  factor  characteristic  of  the 
travel  path.  Determination  of  appropriate  values  for  Q, 
either  in  the  form  of  a Q-depth  profile  or  a characteristic 
T/Q  for  a given  source-receiver  pair  is  an  important  and 
active  research  area. 

A necessary  consequence  of  an  absorption  operator  of 
the  kind  discussed  here  is  the  presence  of  dispersion. 
Theories  to  account  for  this  dispersion  have  been  given  by 
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compared  liy  Savage  and  O'Neill  (1975).  The  difficulty  with 

all  these  theories  is  that  an  arbitrary  cutoff  frequency 

is  required  beyond  which  the  constant  Q assumption  no 

longer  obtains.  Choices  of  this  cutoff  frequency  at  the 

igh  or  low  ends,  along  with  assumptions  about  the  physical 

bases  for  the  cutoff,  account  for  the  differences  between 
the  theories. 

The  three  absorption  operators,  including  attenuation 
and  dispersion,  are  as  follows: 

Lomnitz  (1957); 


exp  {-ituT 


1/2; 


1 ' ttQ  £n  (r1)]  j # exp 
*n  y = 0.577216  ..., 

a - high  frequency  cutoff, 
Futterman  (1962)  ; 


exp  J-iwT 


'-kk**  (yj, 


u0  = low  frequency  cutoff, 


Strick  (1970); 


exp  j- iuT 


i1  • h * bj  (^)]| 

“h  = hiSh  frequency  cutoff. 


oiT 

IQ 


(3.2) 


(3.3) 


(3.4) 


Let  us  compare  a synthetic  seismogram  computed 
with  no  dispersion  to  those  in  which  dispersion  is  included 

using  the  three  formulations  above.  The  cutoff  frequencies 
chosen  are: 
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a = 101 0 , 

« = 10‘5/y,  f3.5) 

0 

to,  = 2000  TT 
h 

The  first  two  are  the  values  suggested  by  Savage  and  O'Neill 
while  the  last  was  selected  by  numerical  experiments  at  S3. 

In  any  event,  the  pulse  shapes  are  quite  insensitive  to  these 
cutoff  values. 

The  seismograms  of  Fig.  3.6  - 3.9  were  computed  using 
the  Strick  formulation.  "Ie  shall  use  that  of  3.6a,  a 10  kT 
explosion  at  a depth  of  0.5  km  in  the  crustal  structure  of 
Table  3.1.  The  epicentral  distance  is  4,000  km  and  T/Q  = 

1.1  was  used.  An  LRSM  instrument  transfer  function  was  ap- 
plied. In  Fig.  3.14,  the  same  seismogram  is  shown  with 
the  Q operator  of  Futterman  and  Lomnitz  being  applied  and 
also  with  no  dispersion  applied. 

Comparing  Figs.  3.14a-d,  wc  first  note  that  the  three 
formulations  of  dispersion  all  have  identical  effects  on  the 
pulse  shape.  In  fact,  when  the  three  pulses  are  overlain, 
they  are  indistinguishable.  Failure  to  include  this  dis- 
persion has  a rather  dramatic  effect  on  the  seismogram,  how- 
ever. 

A peculiar  feature  of  the  dispersion  for  these  linear 
Q operations  is  its  effect  on  the  travel  time.  The  ray 
theoretic  travel  time  is  indicated  on  each  record  by  T. 

The  non-causality  of  the  non-d isper sed  pulse  is  clearly 
evident  in  3.14a  where  the  pulse  begins  before  the  arrival 
time  and  a sharp  pulse  onset  is  not  evident.  Adding  disper- 
sion, a sharp  onset  appears  but  the  pulse  is  moved  several 
seconds  from  the  arrival  time  - later  in  the  high  cutoff  fre- 
quency cases  and  earlier  ir.  the  low  cutoff  frequency  formula- 
tion of  Futterman.  Though  Futterman 's  work  is  often  cited 
by  those  computing  synthetic  seismograms,  the  bothersome 
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20.0 


.799 


r 


6-i sec — l's  i A A 


1.08  L 


< 1 
t 

1 11 


0 5.0  10.0  15.0  20.0 

(b)  Dispersion  formulation  according  to  Lomnitz  (1957) 
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(c)  Dispersion  formulation  according  to  Futterman  (1962) 
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(d)  Dispersion  formulation  according  to  Strick  (1970) 


Figure  3.14.  Comparison  of  synthetic  seismograms  at  4000  krt 
for  attenuation  without  dispersion  and  three  / 
different  formulations  of  the  required  dis- 
persion. The  event  is  the  same  as  that  in 
Fig.  3.6a  and  T/Q  = 1.1.  The  ray  theoretic/ 
travel  time  is  indicated  by  T.  ' 
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forward  time  shift  introduced  by  his  formulation  has  not 
been  mentioned  to  our  knowledge.  However,  this  shift  is 
quite  apparent  in  the  examples  given  in  his  1962  paper. 

The  shifts  to  later  time  in  the  Lomnitz  and  Strick  formula- 
tions are,  of  course,  hardly  more  satisfying. 

The  form  of  the  travel  time  shift  introduced  by  the 
dispersive  operators  (3.2  - 3.4)  can  be  deduced  by  some 
elementary  manipulations.  Let  the  elastic  or  rav  theoretic 
travel  time  be  denoted  T.  Then  the  actual  travel  time  for 
the  three  formulations  is: 

Lomnitz 


(3.6) 


(3.7) 


(3.8) 


For  fixed  values  of  the  cutoff  frequencies  the  travel 
time  shift,  AT,  is  therefore  linearly  proportional  to  T/Q. 

A series  of  synthetic  seismograms  are  constructed  at  T/Q  = 
0.55,  1.1,  2.2,  for  each  of  the  dispersion  formulations  to 
illustrate  this  proportionality.  The  set  at  T/Q  = 11  is 
shown  in  Fig.  3.14.  For  these  cases  the  travel  time  shift 
AT  is  plotted  versus  T/Q  in  Fig.  3.15.  It  should  also  be 
mentioned  that  the  near  identity  of  the  pulse  shapes  for 
all  three  formulations  which  was  apparent  in  Fig.  3.14 
continues  to  hold  for  these  widely  different  T/Q  values. 
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earth,  it  is  also  dependent  on  the  value  ol  the  cutoff  fre- 
quency, a rather  arbitrary  quantity.  Multiplying  the  cutoff 
frequency  by  M,  the  travel  time  shift  changes  by  approximately 
T £n  M/ttQ.  This  is  apparent  by  comparing  the  Strick  and 
Lomnitz  formulations  in  Eq.  (3.6)  and  (3.8)  and  Fig.  3.15. 

The  selected  values  of  the  cutoff  frequencies  are  such  that 

a/y  = 107  co,  . This  accounts  for  the  considerably  larger 
* * 

shift  for  the  Lomnitz  formulation,  yet  the  pulse  shapes  are 
identical.  The  same  kind  of  comparison  can  be  made  using  the 
Futterman  operator.  In  any  case,  the  variation  in  pulse  shape 
among  seismograms  computed  with  cutoff  frequencies  vary- 
ing over  several  orders  of  magnitude  about  the  values  of  Pq . 
(3.5)  is  like  the  trivial  variation  apparent  in  Fig.  3.14 
where  the  three  formulations  were  compared. 

The  time  shift  introduced  by  these  forms  for  including 
dispersion  is  troublesome  because  of  the  ambiguity  it  intro- 
duces in  the  travel  time  and,  therefore,  in  interpretations 
of  travel  time  data.  We  have  three  different  formulations, 
none  of  which  can  be  discarded  on  theoretical  grounds.  They 
have  essentially  identical  effects  on  the  pulse  shape.  Yet 
the  resulting  travel  time  differs  by  as  much  as  tens  of 
seconds  depending  on  T/Q  and  the  cutoff  frequency  chosen. 

That  is,  for  the  same  earth  structure  the  travel  time  varies 
widely  depending  on  the  way  dispersion  is  introduced. 

We  have  posed  an  unresolved  question  which  certainly 
needs  greater  attention.  Fo‘r  the  time  we  have  chosen  to 
use  the  Strick  formulation  tor  routine  applications,  if  only 
because  tne  time  shift  is  minimal  for  this  case. 
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It  is  well  known  that  the  first  few  cycles  of  the 
teleseismic  short  period  P wave  record  of  explosions,  and 
hence  mb>  are  strongly  influenced  by  the  free  surface  re- 
flected phase  pP.  This  effect  has  been  studied  in  previous 
S3  reports  (e.g..  Cherry,  et  ah  (1974b)).  It  has  also  been 
pointed  out  by  a number  of  investigators  that  details  of  the 
near  source  crustal  structure  can  affect  the  shape  of  the 
teleseismic  P wave  record.  A discussion  of  this  effect  was 
given  in  our  last  contract  report  (Bache,  et^  al.,  (1974)). 

A simple  example  will  suffice  to  illustrate  the  kinds  of 
effects  that  may  occur. 

Using  the  theory  of  Appendix  B,  a spherically  symmet- 
ric explosion  source  may  be  buried  in  a layered  model  of 
the  crust.  A series  of  illustrative  synthetic  seismograms 
are  shown  in  Fig.  3.16.  These  are  appropriate  to  an  epi- 
central  distance  of  4, BOO  kin  (one  mantle  arrival  is  computed) 
with  T/Q  = 0.84.  The  source  is  a spherically  symmetric  20 
kT  explosion  in  NTS  tuff  which  is  characterized  by  a source 
spectrum  much  like  #1  of  Fig.  3.2.  The  explosion  is  buried 
at  five  different  burial  depths  in  the  indicated  crustal 
structure  which  is  typical  of  NTS.  Therefore,  the  important 
contributing  phases  may  include  arrivals  from  the  tuff-hard 
rock  interface  at  2.5  km  as  well  as  P and  pi’.  With  each 
seismogram  is  indicated  the  probable  measurement  for  conven- 
tional mb.  The  results  indicate  that  burial  depth  alone  can 
contribute  to  mfe  changes  up  to  0.3  units.  The  coupling  into 
elastic  waves  of  an  explosion  of  constant  yield  will,  of 
course,  be  dependent  on  depth  of  burial  in  a complex  way  as 
the  rock  properties  are  depth  dependent. 

At  the  beginning  of  this  section  in  Fig.  3.1  a set 
of  seismograms  for  five  similar  NTS  events  were  shown.  As 
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source  itself  and  the  crustal  structure  in  the  immediate 
source  vicinity.  A subsequent  discussion  of  the  character 
of  the  one-dimensional  source  characteristics  for  NTS  tuff 
explosions  made  clear  that,  even  over  a broad  yield  range, 
one  must  look  for  higher  order  perturbations  on  the  center 
of  dilatation  to  explain  source  generated  teleseismic  wave 
form  differences  from  event  to  event.  Near  source  crustal 
structure  may,  of  course,  always  play  a part. 

In  order  to  gain  some  insight  into  the  extent  to  which 
fine  details  of  the  crustal  structure  may  be  responsible  for 
the  differences  apparent  in  Fig.  3.1,  an  attempt  to  match 
these  seismograms  using  elementary  models  was  used.  The  re- 
sults of  this  exercise  will  be  presented  here.  Details  of 
the  event  and  receiver  characteristics  are  omitted  for 
security  reasons,  but  these  should  be  inessential  for  our 
purposes . 


In  computing  synthetic  seismograms  to  match  the  data 
of  Fig.  3.1,  the  following  information  was  used: 

1.  Reasonable  estimates  of  the  material  properties 
at  the  shot  point  were  available  allowing  a 
finite  difference  (SKIPPER)  calculation  of  the 
reduced  displacement  potential  (RDP) . Analysis 
of  the  data  indicated  that  one  RDP  scaled  to  the 
proper  yield  would  suffice. 


2.  Along  with  the  event  yield,  the  depth  of  burial, 
average  source-surface  travel  times  and  approxi- 
mate depth  to  the  underlying  Paleozoic  rock 
interface  were  known  for  the  five  events.  Also, 
a detailed  average  velocity-density  profile  for 
the  event  region  to  a depth  of  1 km  was  compiled. 
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3.  Average  crustal  structure's  Tor  NTS  (Hill  ;iiul 
Pakiscr  (1967))  and  the  receiver  region  were 
used.  Since  the  receiver  was  located  on  a 
transparent  region  of  the  crust,  the  receiver 
crustal  structure  had  little  effect. 

4.  The  response  characteristics,  both  amplitude  and 
phase,  of  the  sensing  instrument  were  known. 

5.  For  the  travel  in  the  upper  mantle,  T/Q  = l.l 
was  used.  The  geometric  spreading  was  computed 
with  geometric  ray  theory  (Julian  and  Anderson 
(1968))  and  the  upper  mantle  model  CIT  109  of 
Archambeau  , et^  a 1 . (1969).  Since  geometric 
ray  theory  leads  to  inaccuracies  in  amplitude 
calculations,  the  geometric  spreading  factor  is 
only  approximate.  Therefore,  the  maximum  ampli- 
tudes of  the  observed  and  synthetic  records  for 
event  A were  scaled  by  multiplying  all  synthe- 
tics by  a factor  of  2.1. 

The  results  are  shown  in  Fig.  3.17.  The  agreement  of 
the  shape  of  the  wave  form  for  the  first  few  cycles  of  motion 
is  quite  good. 

The  extended  train  of  higher  amplitude  motion  appearing 
on  the  observed  seismogram  may  be  due  to  later  arriving  mantle 
phases.  While  only  one  mantle  arrival  was  treated  in  our 
calculations  of  four  of  the  five  cases,  an  indication  of  the 
improved  match  which  is  possible  by  including  later  mantle 
phases  is  shown  in  event  D where  a second  arrival  appears  on 
the  synthetic  record. 

Also  indicated  on  each  seismogram  is  a consistently 

measured  mb  using  the  phase  indicated  on  the  event  A on  each 
record . 
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igure  3.17.  A comparison  of  observed  and  synthetic  seismo- 
grams for  five  similar  NTS  events.  A consistent 
mb>  measured  as  indicated,  is  given  for  each 
record . 


74 





U4ilWui  UJ  U J..H  i ... 


I.  -H 


y 

*u 


r 


R-  2646 


3.5 


^iLmIw'amv^1  0F  l,TnillR  0|,n,R  IM-HTIIliHAT  IONS  ON  'll 

MJi'iiilA1'.1-.) ..  symmktr  i c i;x  i,i.o_sij>N_<:i  nisKati  m i-  wav i: 


The  ex  is tenco  of  SH  and  Love  waves  as  well  as  the 
azimuthal  dependence  in  observed  Rayleigh  wave  amplitudes  are 
strong  indicators  that  a simple  center  of  dilatation  model 
of  an  underground  explosion  is  insufficient.  Higher  order 
perturbations  on  the  spherically  symmetric  monopole  source 
are  apparently  present.  A number  of  mechanisms  accounting 
for  these  effects  have  been  postulated  including  tectonic 
release,  earthquake  triggering,  jointing,  spallation,  etc. 

The  theory  r'  Appendix  B was  developed  in  order  to 
provide  an  interface  between  seismic  wave  propagation  methods 
and  the  nonlinear  finite  difference  codes  through  which  many 
of  the  complicated  phenomena  leading  to  a departure  from 
spherical  symmetry  may  best  be  studied.  A few  simple  examples 
which  provide  a qualitative  estimate  of  the  poss.ole  effect 
on  the  teleseismic  short  period  record  will  n0w  be  given  for 

one  important  cause  of  higher  order  source  perturbations, 
tectonic  release. 


Archambeau  (1972)  has  given  a theory  for  the  dynamic 
stress  relaxation  due  to  the  creation  of  a cavity  and  sur- 
rounding zone  of  fractured  material  in  a prestressed  geo- 
logic formation  (tectonic  release).  We  now  consider  an  explo- 
sion plus  tectonic  release  and  examine  the  teleseismic  P 
wave  seismogram. 

The  spherically  symmetric  explosion  source  will  be 
taken  to  be  that  denoted  as  »4  in  Fig.  3.2,  scaled  to  65  kT. 
That  is,  we  take  a rather  weak  coupling  explosion  in  NTS 
tuff.  For  the  tectonic  release  calculation  the  prestress  was 
taken  to  be  120  bars.  The  radius  to  which  shock  induced  frac- 
ture begins  is  Rc  - 30  meters.  This  fracture  then  propagates 
to  a final  rnHiiti;  d - ? nn  


to  a final  radius  of  RQ  = 300  meters.  These  values  fir  the 


tectonic  release  model  compare  to  the  SKIPPER  computed  cavity 
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radius  0!  5-1  meters  and  an  elas.ie  radius  of  motors  , 
physical  parameters  for  the  explosion  and  tectonic  release 
calculations  are  therefore  not  entirely  compatible,  but  this 
is  not  of  much  importance  here. 

A source  consisting  of  an  explosion  plus  tectonic 
release  will  be  composed  of  a center  of  dilatation  plus  a 
ouble  couple.  In  the  notation  of  Appendix  B,  the  nonzero 

an^the  dC  CblffiCientS  Wl11  be  A°>>  for  the  explosion 
and  the  double  couple  or  quadrupole  terms  Ap>  B0). 

m 0,1,2,  j = 1,2, 3, 4;  excited  by  tectonic  release”'  The 
nontero  quadrupole  terms  are  determined  by  the  orientation 
Of  the  prestress  field. 

that  i,LethUS  flrSt  C°nSider  3 hori2°"tally  oriented  prestress 

vert'  th"  PreStr0SS  isi„s  12n  bars,  where  the  3-axis  is 
vertical.  The  geometry  is  depicted  in  Fig.  3.18.  In  this 

case  the  nonzero  multipole  coefficients  arc  for  the 

P wave  excited  by  tectonic  release  and  = AC3J 

“ ' A2/0)  for  the  S wave.  21  22 

In  order  to  visualize  the  elastic  waves  excited  bv 
F1“  ,S°"Ce’SeVeral  Fadiation  pattern  plots  are  given  in 
g.  3.19.  In  the  flgUre  are  several  two-dimensional  sec 
tions  showing  the  1 second  spectral  component  of  the  P and 
SV  wave  displacement  1.0  km  from  the  source.  This  source  is 

°-  C°UrSd:  2 StrOTg  Tadiator  of  Sll  waves,  but  these  are  not’ 
hown.  Three  of  the  plots  show  the  vertically  emitted  radia- 

5 Ct6d  azi"uths-  shall  subsequently  present 
synthetic  seismograms  computed  at  4,000  km.  At  this  distance 
e earth  model  used  gives  a takeoff  angle  of  14.7°  for  the 
direct  P wave  and  7.8°  for  the  free  surface  reflected  ,P 
Phase.  The  takeoff  angles  for  the  three  main  phases,  P,  PP 

earthrmodel1Caded  ^ ^ aZi"Uthal  Pl°tS-  F°r  a"^  ™«onabie 
model  and  over  the  entire  teleseismic  distance  range 
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Figure  3.18 . 


Coordinate  system  for  defining  the  nrestr^ 

Therfirstdno‘S°Ur^0"tu"reCeiver  azin,uth  (0) 
me  r irst . motion  m the  four  quadrants  fi,L 

an  explosion  m a medium  with  prestress 

and  a ^dicatcd*  ^finUions  of  , 

and  a23  are  consistent  with  that  of  0 

1 2 * 


77 


C J 


revise3 radiation  patterns  for  an  explosion  plus  tectonic 
release.  Depicted  are  the  1 Hz  spectral  components  of  the  P and  SV  wave 

J,w.ta  fromuth®  source.  The  contribution  of  the  spherically 
!^ftric  exPl°slon  to  the  P wave  is  indicated  separately  by  the  dashed  Y 

vUTVC  • 
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Figure  3.19  (Continued) 
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the  takeoCf  angles  of  interest  arc  not  far  Jiff,. rent  from 
these.  Some  idea  of  the  azimuthal  pattern  is  giver.  .,y  the 

fourth  plot  in  which  is  depicted  the  P and  SV  waves  omitted 

at  the  typical  teleseismic  takeoff  angle, G = 15°. 

It  is  clear  that  the  orientation  of  the  prestress  will 
strongly  influence  the  ratio  of  the  P and  S waves  directed 
into  the  teleseismic  field.  It  is  useful  to  look  at  the 
orientation  at  which  the  greatest  enhancement  of  S waves 

will  occur;  that  is,  the  orientation  at  which  the  S wave 

lobe  is  directed  vertically.  Such  an  orientation  is  given 
by  rotating  the  prestress  field  o;  Fig.  3.18  by  45°  about 
the  x2  axis,  then  by  45°  about  the  axis.  The  result  is 
a prestress  field  that,  in  earthquake  parlance,  has  strike 
dip  and  plunge  of  90°,  45°,  45°.  Using  the  theory  worked  ’ 
out  by  Minster  (1974)  for  rotating  the  multipole  coefficients 

WP  f i nrl  t-  Vs  rs  a.  . ~ i - r a \ 


o - "‘v-*  i.  x V L.U  C I 

we  fmd  the  new  nonzero  quadrupole  coefficients  are  A 
n ? c * (4)  , „ . ms 


(4) 


°'n)A22  ^ r"^'5  A 2 1 for  the  P wave  and  (w)°  = 

3R(-iJr„i'l  = . n (3)  r..\  1 ^.(2).  . __fn  20  ri^J 


(“0 


3BCl)(m)  = - B_C3)(o))  =*-  6A.C2)(u>)  = - 2B(1)  (u,)  = 6B^  (w) 


(3) 


21  21  22''  — . . v~  J (w  j 

for  the  S wave.  For  this  orientation  the  one  second’spectra 
component  of  the  S wave  displacement  is  about  4.0  times  as 
large  as  the  corresponding  P wave  displacement  along  the 
teleseismic  takeoff  angles. 

We  have  now  constructed  two  sources  consisting  of  an 
explosion  in  tuff  plus  tectonic  release.  The  tectonic  re- 
lease component  is  rather  large,  perhaps  providing  an  upper 
bound  to  what  may  be  encountered.  In  one  case,  the  prestres* 
field  is  horizontal,  while  in  the  second  it  is  rotated  to 
maximize  the  sP  phase  on  the  teleseismic  record.  Let  us 
characterize  the  two  sources  by  the  ratio  of  the  displacement 
along  the  direction  marked  P to  the  direction  marked  sP  in 
the  takeoff  angle  radiation  patterns  of  Fig.  3.14.  This 
ratio  will  be  denoted  P/SV. 
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l,;'vinK  1 onst  rue  t oil  t ho  so  sources,  let  u-.  now  Mn.ly 
the  result  ini’,  tcieseismic  short  period  I1  wave  record.  l.et 
the  explosion  be  buried  at  a depth  of  400  meters  in  a typi- 
cal NTS  crustal  structure.  A series  of  seismograms  at  an 
epicentral  distance  of  4,000  km  are  shown  in  Figs.  3.19  - 
3.20  where  T/Q  = l.l  was  used. 


In  Fig.  3.20  two  crustal  structures  are  used,  one 
with  the  tuff-hard  rock  interface  at  2 . 5 km  and  one  with  it 
at.  0.7  km.  Seismograms  for  both  sources,  that  with  the  pre 
stress  horizontal  (P/SV  = 2.0)  and  that  for  which  the  pre- 
stress field  was  rotated  to  optimize  the  sP  phase  (P/SV  = 

O. 25),  were  computed.  The  top  seismogram  shows  the  effect 
of  ignoring  reverberations  in  the  crust;  that  is,  only  the 

P,  pP  and  sP  phases  are  included.  For  the  P/SV  = 2 source 
the  sP  phase  is  too  small  to  have  any  appreciable  effect  on 
the  record.  In  the  second  synthetic  seismogram  crustal 
reverberations  are  included,  having  a major  effect  on  the 
pulse  shape  and  mb< 


The  second  and  third  records  of  Fig.  3.20  give  a 
direct  comparison  of  the  two  sources.  Only  the  orientation 
of  the  prestress  field  differs  between  the  two  but  the  change 
in  mb  of  nearly  0.2  units  is  significant.  Finally,  the  last 
record  indicates  the  possible  importance  of  the  tuff-hard 
rock  interface.  Note  the  deep  second  trough,  a peculiar 

characteristic  of  several  of  the  observed  records  shown  in 
Fig.  3.1. 


In  Fig.  3.21  we  compare  the  pure  explosion  to  the 
explosion  plus  tectonic  release.  The  tectonic  release  com- 
ponent is  expected  to  represent  an  upper  bound  for  the  con- 
tribution to  the  tele seismic  short  period  P wave  record. 
Three  slightly  different  crustal  models  are  compared;  model 
3 of  the  previous  figure,  the  same  model  with  the  tuff -hard 
rock  interface  at  1 km,  and  a model  in  which  the  velocity 
transition  is  less  abrupt.  The  difference  between  the  three 
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P Wave  Velocity  (kni/sec) 


Structure  A 


Structure  B 


Figure  3.20 


Synthetic  seismograms  at  an  cpiccntial  distance 
of  4,000  km  for  an  explosion  plus  tectonic  re- 
lease in  NTS  tuff.  Attenuation  was  characterized 
by  T/Q  = 1.1  and  the  explosion  depth  of  burial 
was  0.4  km.  The  indicated  m^  values  were  com- 
puted in  the  conventional  way  from  the  maximum 


rw*  ~rrr. 


? 


«ii 


R-2646 


cases  is  not  dramatic,  but  indicates  that  the  presence  of  a 
sliaip  velocity  discontinuity  is  more  important  than  its 
exact  location. 


Finally,  these  few  examples  seem  sufficient  to  con- 
clude that  only  under  optimal  circumstances  can  tectonic 

release  have  a noticeable  effect  on  mb  and  the  teleseismic 
P wave  record 
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lv-  VAR  I A I'.  1,1,  I'RI'lllllMCV  MAUN  I Till)!  H I M 'I'  I M I N AMI 
4 . 1 INTRODUCTION 

During  the  past  several  months  we  have  been  actively 
engaged  in  the  development  and  testing  of  techniques  for 
improved  discrimination  between  earthquakes  and  both  single 
and  multiple  underground  explosions.  As  reported  previously 
(Savino  and  Archambeau,  1974;  Bache,  et  aL,  1974),  this 
program  has  resulted  in  the  formulation  of  an  extremely 
promising  and  very  effective  variable  frequency  magnitude 
(VFM)  discriminant  which  exploits  spectral  differences 
between  earthquakes  and  explosions.  To  date  our  attention 
has  focussed  on  the  application  of  the  VFM  technique  to 
short  period  body  waves  from  events  at  teleseismic  distances 
from  recording  sites.  As  will  be  shown,  obviating  the  require 
ment  for  recording  small  explosions  results  in  an  extended 
magnitude  range  over  which  the  VFM  discriminant  can  be  ap- 
plied. Complete  separation  of  earthquake  and  explosion 
populations  is  achieved  down  to  event  magnitudes,  m,  , in  the 
raige  4.0  to  4.  5.  b 

Following  a brief  description  of  the  VFM  technique, 
the  results  of  its  application  to  a large  population  of 
shallow  and  deep  focus  events  recorded  at  LASA  and  the 
original  Oyer  subarray  in  Norway  will  be  described.  Finally, 
the  VFM  technique  will  be  tested  for  its  effectiveness  as 
a discriminant  on  a simulated  multiple  explosion  similar  to 
the  scenario  proposed  by  Kolar  and  Pruvost  (1975). 

4.2  TECHNIQUE 

The  variable  frequency  magnitude  discriminant,  original- 
ly proposed  by  Archambeau,  et  al_^  (1974),  is  designed  to 
exploit  spectral  differences  between  earthquakes  and  under- 
ground explosions.  Essentially,  the  VFM  technique  consists 


1 


85 


4 


~l“r: 


R - 2 6 4 6 


of  a comnr-ison  ol  a magnitude'  men  mi  remon  1 iii(  ||  i ,, 
relatively  low  frequency  (e.g.,  f^  « (l.lf,  II  ) to  a mngni 
tude  measurement  n^ff  ) at  a higher  frequency  fe.j ...  f . 
2.25  Hz).  The  body  wave  magnitude  Eb  is  defined  as  the  log 
of  the  amplitude  of  the  output  from  a narrow-band  (high  Q) 
phaseless  filter,  centered  at  frequency  f £ , plus  a distance 
correction  factor  F.  That  is 


where  the  center  period  T,  = 1/f.  and  ^ is  the  maximum 
amplitude  of  the  envelope  of  the  filter  output . The  filters 
employed  are  digitally  constructed  in  either  the  time  do- 
main (recursive)  and/or  frequency  domain.  The  recursive 
filter  is  made  phascless  by  filtering  an  original  time 
series  in  both  the  forward  and  backward  time  sense  and  sum- 
ming the  two  results. 

Examples  of  narrow-band  recursive  filter  outputs  for 
a presumed  Eurasian  explosion  and  a shallow  earthquake  re- 
corded by  the  Oyer  subarray  in  Norway  are  shown  in  Fig.  4.]. 
The  top  traces  on  the  right-  and  left-hand  sides  of  this 
figure  correspond  to  the  unfiltered  best  beam  recordings  of 
the  P-wave  trains,  preceded  by  about  30  seconds  of  background 
noise,  from  tne  earthquake  and  presumed  explosion,  respec- 
tively. Note  the  enhancement  of  the  high  frequency  (f  =6.0 
Hz)  filter  output  versus  the  low  frequency  (f  =0.3  Hz)  fil- 
ter output  for  the  presumed  explosion  signal  as  compared  to 
the  filter  outputs  for  the  earthquake  signal.  The  m (f  ) 
estimates  for  these  events  would  be  based  on  the  maximum 
amplitudes  of  the  envelopes  of  the  filter  outputs  at  the 
corredponding  f^s.  As  an  example,  the  ^ (6.0  Hz)  estimate 
for  the  presumed  explosion  would  be  based  on  the  amplitude 
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Figure  4.1.  Examples  of  narrow  band  filter  outputs  at 
center  frequencies  CO. 3 Hz,  1.0  Hz  and  6.0  Hz)  for  an 
sion  (left-hand  side)  and  a shallow  earthquake  (right- 
side)  . 
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30  seconds. 

4.3  SHALLOW  EVENTS  RECORDED  AT  LASA 

The  data  base  employed  in  the  first  test  of  the  VFM 
discriminant  consists  of  LASA  short-period  recordings  of 
P-wave  trains  from  34  presumed  explosions  and  156  earthquakes 
(Lacoss,  1969).  Thirty  of  the  presumed  explosions  originated 
within  the  mainland  USSR,  two  in  Novaya  Zemlya,  one  (Long- 
shot)  in  the  Aleutians  and  one  in  the  Sahara  Desert.  The 
earthquakes  are  distributed  along  the  Alpide  seismic  belt, 
the  Kuril -Kamchatka  arc  and  the  Arctic  Ocean.  The  epicentral 
distances  of  these  events  range  from  45°  to  about  100°  with 
more  than  half  of  the  events  between  65°  and  85°  from  LASA. 

Figure  4.2  is  a plot  of  spectral  magnitude  estimates, 
m^f  ),  at  a low  frequency  (fc  * 0.45  Hz)  versus  a high 
frequency  (f£  = 2.25  Hz)  for  all  34  presumed  explosions  and 
those  earthquakes  in  the  data  base  that  were  repo  ed  in  the 
Monthly  Listings  of  Fvents,  published  by  the  USOS  (formerly 
NOAA) , as  having  shallow  focal  depths,  h *"  "O  km.  This  figure 
clearly  demonstrates  the  enriched  high  frequency  content  of 
the  explosion  body  wave  spectra  as  compared  to  the  earthquake 
spectra.  For  instance,  for  a given  value  of  in^  (0.45  Hz) 
the  explosions  exhibit  in^  (2.25  Hz)  values  that  are  typically 
0.6  to  1.0  unit  larger  than  the  m^  (2.25  Hz)  values  for 
earthquakes . 


The  high  degree  of  discrimiantion  of  earthquakes  from 
* explosions  evident  in  Fig.  4.2  is  especially  significant  in 

view  of  the  non-regionalization  of  the  event  population. 

The  variety  of  tectonic  settings  of  this  event  population 
ranges  from  relatively  stable  shield  regions  to  seismically 
p active  oceanic  arc  systems.  An  indication  that  discrimination 

could  be  further  enhanced  by  regionalizing  the  event  popula- 
tion comes  from  the  fact  that  for  in,  (2.25  Hz)  > 4.0  the  two 


88 


i. 


mbC2.25HzD 


Figure  4 
and  2.25 
occurred 


2.  Spectral  magnitudes,  m^,  computed  at  0.45  Hz 
Hz.  The  presumed  explosions  numbered  55  and  138 
at  Novaya  Zcmlya. 
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presumed  "xpl  os  i nn<  (#35  and  M3H)  plotting  closest  |„  the 
oa  rt hquake  populat  ion  occurred  at  \ovay.i  leiiilya.  Several  of 
the  earthquakes  that  plot  closest  to  the  explosion  population 
occurred  along  the  Kuri  1 e - Kamchatka  arc  or  at  locations  far 
removed  from  the  explosion  epicenters. 

The  apparent  bending  of  the  explosion  population  into 
the  earthquake  population  at  m^  (2.25  Hz)  < 3.5  in  Fig.  4.2 
is  mainly  a result  of  microseismic  noise  inflating  the  low 
frequency  (fc  = 0.45  Hz)  magnitude  estimates  for  the  rela- 
tively small  (signal -to-noi se  wise)  explosions.  In  order  to 
see  the  effects  of  noise  contamination  more  clearly,  magni- 
tude estimates  based  on  a series  of  1 ow- frequency  (0.3  to 
0.6  Hz)  narrow  band  filters  versus  the  f = 2.25  Hz  high  fre- 
quency filter  for  a subset  of  the  event  population  in  Fig. 

4.2  are  shown  in  Fig.  4.3a-d.  Note  how  the  trends  of  the 
explosion  population  (x's)  in  the  case  of  the  low  frequency 
(fc  = and  0.4  Hz)  filters,  which  are  closest  to  the  main 

concentration  of  noise  power  at  LASA  between  about  0.1  to 
0.3  Hz  (Lacoss  and  Toksoz.  1967),  exhibit  rather  abrupt  bends 
near  m^  (2.25  Hz)  = 5.0  and  flatten  out  at  mb  (2.25  Hz)  < 5.0. 
Examination  of  Fig.  4.3a-d  indicates  that  the  prevailing  level 
and  spectral  distribution  of  microseismic  noise  at  LASA  sets 
gradually  decreasing  lower  limits  on  the  low-frequency  magni- 
tude estimates  as  the  narrow  band  filter  frequency  increases 
from  0.3  Hz  to  0.6  Hz,  thereby  moving  away  from  the  main 
noise  band.  In  the  case  of  Fig.  4.3d,  for  mb  (0.6  Hz)  versus 
mb  (2.25  Hz),  the  trend  of  the  explosion  population  is  more 
nearly  linear  and  parallel  to  the  earthquake  population  over 
the  entire  magnitude  range  of  the  data.  Noise  contamination 
of  the  earthquake  magnitude  estimates  in  Fig.  4.3a-d  is 
obviously  not  as  serious  a problem  since  the  earthquake  sig- 
nals for  a gi\en  value  of  mb  (2.25  Hz)  are  richer  in  low 
frequencies  and  shou  only  a slight  tendency  to  decrease  in 
value  as  the  filter  frequency  incr  uses  from  0.3  I-  to  0.6  Hz. 
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Fig.  4.3c 
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Fig.  4.3d 


thA^ACA*3*  Spectral  magnitude  estimates  for  a subset  of 
S!  LASA  event  population  showing  the  effect  of  varving 
the  tc  of  the  low  frequency  magnitude  estimates. 
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» the  high  frequency  end  of  the  seismic  bandpass  „ 
the  interplay  of  me  inly  three  factors  that  Jimit  the 
pi icabi 1 it y of  the  VFM  discriminant,  as  well  as  other  short 
period  discriminants,  to  small  magnitude  teleseismic  events 
recorded  at  LASA.  These  three  factors  are  source  spectrum' 

:r  ;:“C  “!enuati0n  and  "^-frequency  background  noise  ’ 
“ f D\fl"ltlVe  inf°™ation  the  source  parameters  of 
earthquakes  included  in  this  data  base  are  not  available 
wever,  the  NOAA  magnitudes  of  all  the  reported  events  are 
greater  than  mb  =■  4.7.  Thus,  we  will  as  sume  that  the  source 

that tfh  C°rner  fTeqUenCieS  than  2.0  Ha  and  further 

(w-i  to6  -mT  SP6Ctra  3re  int°  thC  hiKh  fre‘>uenc>’  roll-off 
( t0  “ ’ f°r  the  hi8»  frequency  range  of  interest  here. 

anelastic  ^ kn°K":  3 flXCd  — ce-receiver  distance 

• . GnUatl0n  glVGS  rise  t0  an  exponential  decay  with 

increasing  frequency  of  the  amplitudes  of  seismic  waves 

us  te  amplitudes  of  P-waves  of  frequency  higher  thanks 

and  d t many  th°  6VentS  ^ the  Particular  magnitude 

: afn::granee  ofinterest  herc  ^ 

tne  nigh  frequency  noise  level  at  LASA. 

The  combined  effect  of  these  different  factors  on  the 
magnitude  estimates,  ff  can  be  seen  in  Fig  4 4,.f  ° \ 

Of  these  figures  the  low  frequency  magnitude  e.im^es  are 

Starting  V **”  ^ 

g nn  tne  tc  1.75  Hz  estimates  in  Fig.  4. 4a  we 
see  that  the  explosion  and  earthquake  populations  are’not 
para  e as  much  as  m Fig.  4.2.  This  is  undoubtedly  due  to 
relatively  close  spacing  of  the  high  and  low  center  fre- 
encies,  0.S  „t  and  1.7S  Ha.  As  the  center  frequency  of  he 
hrgh  frequency  magnitude  estimates  increase,  fFigs  4 4b-dl 
e separation  Of  the  event  population  increaL*  p!^  ’ 

as  the  center  frequency  is  further  increased  fFig:  . 4.4e  and 
f)  the  event  populations  begin  to  merge.  The  majority  of 
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Fig.  4.4a 
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Fig.  4.4b 
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Fig.  4 . 4d 


Figure  4.4.  Spectral  magnitude  estimates  for  a subset  of 
the  LASA  event  population  showing  the  effect  of  varying 
the  fc 's  of  the  high  frequency  magnitude  estimates. 
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ml»(,ci  “ 4,0  ,,t  fc  “ 2-r*  llz  wliile  a comparable  noise  limit  tor 
the  explosions  is  not  reached  until  = 3.0  llz.  Thus,  while 
optimum  discrimination  would  result  for  spectral  magnitudes 
computed  at  frequencies  (fj  separated  by  a decade  or  more, 
noise  and  anelastic  attenuation  of  the  earth  combine  to  con- 
strain the  useable  bandwidth  of  frequencies  in  the  case  of 
the  LASA  data  set  to  approximately  one-half  a decade. 

4-4  DEEP  EVENTS  RECORDED  AT  LASA 

Twenty  events  in  the  depth  range  80  km  to  580  km  are 
included  in  the  LASA  data  base.  Magnitude  estimates  for 
these  events  based  on  narrow  band  filter  outputs  at  the  same 
center  frequencies  as  in  Fig.  4.2  are  plotted  in  Fig.  4.5 
and  compared  with  the  shallow  earthquake  and  explosion  popula- 
tions which  are  contoured  by  the  dashed  and  solid  lines, 
respectively.  All  of  the  events  with  focal  depths  h > 300 
km  occurred  beneath  the  Japan  and/or  the  Kurile-Kamchatka 
arcs.  The  epicentral  locations  of  the  intermediate  depth 
events  range  from  the  Aegean  Sea,  Rumania,  Hindu  Kush  regions 
to  the  Japan  and  Kurile-Kamchatka  arcs. 

Deep  earthquakes  often  fail  to  separate  from  explosion 
populations  when  examined  with  any  of  the  discriminants, 
short-and/or  long-period,  proposed  ro  date.  As  can  be  seen 
in  Fig.  4.5,  such  is  the  case  for  several  of  the  deep  events 
examined  with  the  VFM  technique.  Suggestions  as  to  the  rea- 
son for  the  lack  of  discrimination  of  many  deep  events  range 
from  impulsive  source-time  functions,  small  source-dimensions, 
and  high-Q  propagation  path  to  the  receiver.  While  the  be- 
havior of  the  deep  events  in  Fig.  4.5  cannot  be  explained  at 
this  time,  it  should  be  noted  that  the  majority  of  intermediate 
depth  events  lie  within  the  bounds  of  the  shallow  earthquake 
population. 
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4.5  Si  [ALLOW  I.VFNTS  RFCOROLD  AT  NORWAY 

A subset  of  the  event  population  recorded  at  I ASA 
was  also  recorded  by  the  limited  Oyer  array  in  Norway.  The 
spectral  magnitude  discriminant  was  applied  to  this  smaller 
population  of  events  and  the  results  are  presented  in  Figs. 

4.6  and  4.7.  In  Fig.  4.6,  spectral  magnitudes,  m^,  at  0.6 

Hz  and  5.0  Hz  are  givjn  for  shallow  earthquakes  (open  circles) 
and  presumed  explosions  (closed  circles).  Analogous  to  the 
results  at  LASA  in  Fig.  4.7,  there  is  discrimination  of 
events  over  most  of  the  magnitude  range  of  the  events  except 
for  the  smallest  magnitude  explosion. 

The  arrows  attached  to  the  explosion  points  indicate 
that  the  low  frequency  magnitude,  inb  (0.6  Hz),  is  contaminated 
by  noise.  In  Fig.  4.7  the  output  amplitudes  from  the  narrow- 
band  filters  are  corrected  for  noise  and  replotted.  Comparing 
Figs.  4.6  and  4.7  it  becomes  clear  that  the  cause  of  the  bend 
in  the  explosion  population  toward  the  earthquake  population 
at  small  magnitudes  is  the  inflation  of  the  low  frequency 
explosion  magnitudes  by  noise.  Correction  for  this  noise  re- 
sults in  complete  separation  of  the  earthquake  and  explosion 
populat ions . 

The  effects  of  noise  in  different  frequency  bands  can 
be  seen  in  Fig.  4.8.  In  this  figure  the  unfiltered  time 
series  corresponding  to  the  smallest  explosion  recorded  at 
Norway  is  plotted  in  the  top-left  frame.  The  remaining 
frames  slow  the  effect  on  the  time  series  of  the  application 
of  five  narrow  band  filters  of  increasing  center  frequency 
(0.3  Hz  to  6.0  Hz).  The  gradual  emergence  of  the  signal  at 
fc  > 1*0  Hz  is  quite  striking.  The  persistance  of  a high 
signal -to-noise  ratio  for  this  very  small  event  probably  re- 
sults from  the  high-Q  nature  of  the  propagation  path. 
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Figure  4.7.  Filter  amplitudes  (maximum)  with  noise  correc- 
tions for  the  same  event  population  plotted  in  Fig  4 6 
Note  the  enhanced  separation  of  the  earthquake  and ’explosion 
populations. 
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4.6  nn  r i-vi'nts  RK:oi'’i)i  h_at  Norway 

hi even  events  in  the  Norway  data  base  were  assigned 
focal  depths  > 70  km.  Amplitudes  of  narrow-band  filter  out- 
puts with  the  same  f c ' s used  in  Figs.  4.7  and  4.8  were  com- 
puted for  these  eleven  events  and  are  plotted  in  Fig.  4.9. 

The  frequency  dependent  amplitudes  were  corrected  for  noise 
based  on  measurements  taken  before  the  onset  of  each  event 
P-wave  signal.  The  bounds  of  the  shallow  earthquake  and  ex- 
plosion populations  plotted  ir  Figs.  4.7  and  4.8  are  in- 
dicated in  Fig.  4.9  by  the  broken  and  dashed  curves,  respec- 
t ively . 

As  was  the  case  with  the  deep  events  in  the  LASA  data 
set  shown  previously  in  Fig.  4.5,  the  deep  events  recorded 
at  Norway  (Fig-  4.9)  exhibit  considerable  scatter  without 
any  obvious  pattern  of  behavior.  While  several  of  these 
deep  earthquakes  fail  to  discriminate  in  this  figure  they  are 
really  not  as  troublesome  as  they  might  seem  since  these 
events  can  often  be  identified  as  naturally  occurring  earth- 
quakes on  the  basis  of  other  information  (e.g.,  hypocentral 
location  using  P,  pP  and/or  sP  travel  times). 

A rather  dramatic  example  of  the  short-period  spectral 
similarities  between  an  explosion  and  a deen  earthquake 
is  shown  in  Fig.  4.10.  Approximately  80  seconds  of  an  un- 
filtered explosion  time  series  and  the  output  amplitudes  of 
three  narrow  band  filters  (f.  = 0.3  Hz,  1.0  Hz  and  6.0  Hz) 
applied  to  the  time  series  are  shown  on  the  left-hand  side  of 
this  figure.  A similar  sequence  of  pictures  is  shown  on  the 
right-hand  side  for  a deep  (h  = 115  km)  earthquake.  Note 
how  the  signal-to-noise  ratio  of  the  narrow  band  filter  out- 
puts for  both  events  increases  as  the  center  frequencies  of 
the  filters  increase  from  03.  Hz  to  6 Hz.  Comparison  of  the 
spectral  behavior  of  the  deep  earthquake  in  this  figure  with 
the  shallov;  earthquake  (h  = 35  km)  shown  in  Fig.  4.1  clearly 
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Figure  4.9.  Amplitudes  of  narrow  band  filter  outputs  at  two 
different  frequencies  for  deep  earthquakes  recorded  it  Norway. 
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Figure  4.10.  Narrow  band  filter 
center  frequencies  ( f c = 0.3  Hz, 
explosion  (left-hand  side)  and  a 
side) . 
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4 . 7 MULTIPLE  EXPLOSION  SCFNARIO 

The  spectral  magnitude  technique  is  especially  suited 
ror  identification  and  discrimination  of  multiple  explosion 
sequences  that  are  designed  to  appear  earthquake -1 ike  in 
terms  of  conventional  (Ms-mbf  depth  of  focus,  complexity, 
first  motion)  discriminants.  A multiple  event  scenario,  some- 
what similar  to  one  proposed  by  Kolar  and  Pruvost  (1975),  was 
devised  by  superposing  eight  scaled  seismograms  of  a presumed 
Kazakh  explosion  recorded  at  LASA.  This  explosion  signature  is 
shown  on  the  bottom-center  in  Fig.  4.11. 


The  array  of  explosions  and  their  relative  yields 
were  designed  to  produce  earthquake- 1 ike  seismograms  over 
a wide  range  of  azimuths.  The  particular  array  configuration 
(spacing  and  firing  order)  is  indicated  in  the  center  of 
Fig.  4.11.  Each  seismogram  comprising  the  multiple  event 
was  delayed  in  time  relative  to  the  first,  and  scaled  in 
amplitude.  The  largest  explosion  in  the  scenario  is  the  sixth 
event  and  is  scaled  to  give  the  same  teleseismic  ground  motion 
as  the  primary  signal.  The  amplitude  scaling  for  all  eight 
explosions  in  th«  scenario  was  1,  3.1,  5,  10,  10,  20,  15  and 
12.5. 


The  composite  seismograms  resulting  from  the  scenario 
are  shown  in  Fig.  4.11  at  five  different  azimuths  (1-5)  with 
respect  to  the  shot  array.  The  first  point  to  be  noted  about 
these  composite  seismograms  is  that  with  the  addition  of 
noise  to  the  beginning  of  each  seismogram  the  first  motions 
at  azimuths  1 , 3 and  5 would  most  probably  be  picked  as  rare- 
factions. Secondly,  the  complexity  of  each  composite  signal 
has  been  greatly  increased  over  that  of  the  primary  signal. 
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'.con.irio  .1  p|>(';i  r r:i  i t lupi.i  kc  - J i ko  on  ;m  Ms  m h.r  i..  rMS, 
analyst  making  amplitude  measurements  to  be  used  for  m de- 
terminations in  the  conventional  manner  (maximum  amplitude 
within  the  first  3 or  4 cycles)  would  undoubtedly  pick  ampli- 
tudes corresponding  to  the  earlier  smaller  explosions  in  the 
sequence.  On  the  other  hand  the  Mg  measurements  would  be 
mainly  based  on  the  superposed  surface  waves  from  the  three 
large  yield  explosions  occurring  late  in  the  sequence.  The 
iiet  result  would  be  a reduced  mb  and  an  enhanced  M result- 
ing m the  scenario  moving  into  the  earthquake  population 
on  an  Ms-mb  plot. 

Application  of  the  spectra]  magnitude  technique,  how- 
ever, results  in  complete  discrimination  of  the  multiple 
explosion  scenario.  In  Fig.  1.12  the  spectra]  magnitude  es- 
timates, computed  for  the  five  azimuths  indicated  in  Fig. 

4.11  are  indicated  by  the  x’s.  A most  significant  result  in 
Fig.  4.12  is  that  the  spectral  magnitude  estimates  of  the 
multiple  explosion  sequence  cluster  around  the  estimates  of 
the  primary  signal  used  in  the  construction  of  the  scenario 
(the  closed  circle  immediately  to  the  right  of  x-1). 

This  means  that  the  7FM  technique  has  based  the  magnitude 
estimates  on  the  largest  amplitude  arrival  in  the  wavetrain, 
corresponding  to  the  largest  yield  explosion  in  the  sequence. 
This  is  a very  important  result  for  yield  determination  of 
both  single  and  multiple  explosions. 


4.8  SUMMARY 

A variable  frequency  magnitude  technique  designed  to 
exploit  spectral  differences  between  earthquakes  and  explo- 
sions was  developed  and  applied  to  a large  population  of 
Eurasian  events  recorded  at  LASA  and  a limited  array  (Oyer) 
m Norway.  The  magnitude  estimates  are  based  on  the  output 
amplitudes  of  variable  frequency  narrow  band  phaseless  filters 
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* igure  4.12.  Spectral  magnitudes  for  same  event  population 
and  filter  parameters  as  in  Fig.  4.2  with  estimates  for 
multiple  explosions  (x's)  showing  complete  discrimination. 
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applied  *o  short  period  recordings  of  body  waves.  The  most 
significant,  results  pertaining  to  event  identification,  dis- 
crimination and  yield  determination  obtained  with  this  dis- 
criminant are  the  following: 

1.  Complete  and  positive  discrimination  of  earth- 
quakes and  explosions  down  to  explosion  magni- 
tudes of  LASA  mb  = 4.2. 

2.  Multi-azimuth  discrimination  of  a multiple 
explosion  scenario  with  the  correct  identifi- 
cation of  the  largest  yield  explosion  in  che 
scenario . 

3.  Simultaneous  discrimination  at  a single  record- 
ing site  of  events,  both  earthquakes  and  explo- 
sions, distributed  over  a wide  geographical  re- 
gion. 
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V.  MULTIPLE  DISCRIMINANT  TI-SM 

5. 1 INTRODUCTION 

During  the  past  10  years  or  more,  several  different 
techniques  for  discriminating  between  earthquakes  and  under- 
ground explosions  have  been  proposed.  These  techniques, 
which  have  met  with  varying  degrees  of  success,  depend  upon 
differences  in  the  time-  and/or  frequency-domain  between 
earthquake  and  explosion  signals.  In  this  section  of  the 
report  we  will  compare  the  effectiveness  of  the  VFM  dis- 
criminant to  the  following  proposed  discriminants;  complexity, 
spectral  ratio,  higher  moments  of  frequency  and  m^-M  . A 
most  important  aspect  of  this  comparison  is  that  all  the 
discriminants  are  applied  to  the  same  data  base:  the  large 

Eurasian  population  of  events  recorded  at  LASA. 

5.2  COMPLEXITY 

In  the  early  1960's  the  research  group  at  Blacknest 
proposed  that  a measure  of  the  complexity  of  the  P-wave 
trains  from  earthquakes  and  underground  explosions  could 
possibly  serve  as  a technique  for  identifying  these  two 
types  of  events  (Thirlaway,  1963;  Carpenter,  1963  and  1964). 
They  noted  that  telesei smically  recorded  P-wave  signals  from 
underground  explosions  were  often  simple,  consisting  of  one 
or  two  cycles  of  relatively  large  amplitude  followed  by  a 
tail  or  coda  of  much  lower  amplitude.  P-wave  signals  from 
earthquakes,  however,  often  consisted  of  a series  of  arrivals 
of  comparable  amplitude  spread  over  several  tens  of  seconds. 
The  original  definition  of  complexity  (Carpenter,  1964)  in- 
volved comparison  of  the  inverse  ratio  of  the  signal  energy 
in  the  first  five  seconds  of  the  P-wave  train  to  the  signal 
energy  in  the  subsequent  30  seconds.  This  definition  was 
applied  to  the  cross-correlation  function  between  two  steered 
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In  the  late  1960's,  members  of  the  Lincoln  Laboratory 
research  group  (Kelly,  1968;  Lacoss,  1969)  proposed  a modi- 
fied method  of  computing  signal  complexity  and  applied  it 
to  Eurasian  events  recorded  at  LASA.  In  particular,  Lacoss 
made  linear  complexity  measurements  on  best  LASA  beams  for 
bandpass  filtered  (0.6  Hz  to  2.0  Hz)  P-wave  signals  using 
the  following  definition  of  complexity: 


C 


L 


5 


x(t)  | dt/ 


x (t)  [ dt 


(5.1) 


where  x(t)  is  the  time  series  corresponding  to  the  best 
filtered  beam,  the  time-origin  is  chosen  at  the  event  onset, 
and  in  practice  the  integrals  are  replaced  by  summations  in 
the  case  of  digitized  data  samples. 

Using  Fq.  (5.1),  values  of  complexity  were  computed 
for  the  same  population  of  shallow  earthquakes  and  explosions 
plotted  in  Fig.  4.2.  These  complexity  values  are  plotted  in 
Fig.  5.1  versus  mfe  (LASA).  While  many  of  the  earthquakes 
separate  by  nearly  an  order  of  magnitude  ^rom  the  explosion 
population,  there  is  significant  overlap  of  the  two  popula- 
tions over  the  entire  magnitude  range.  Thus,  in  terms  of 
complete  separation  of  event  populations  it  is  quite  obvious 
that  the  V FM  discriminant  (Fig.  4.2)  i&  superior  to  complexity. 

The  behavior  of  complexity  in  the  case  of  very  small 
magnitude  events  is  controlled  by  background  noise.  As  the 
magnitudes,  and  hence  the  signal-to-noise  ratios,  of  events 
decrease,  the  values  of  complexity  for  both  earthquakes  and 
explosions  (Fig.  5.1)  approach  a value  of  CL  = 6.  This  is 
to  be  expected  since  in  the  limiting  case  of  pure  noise 
(assumed  to  be  stationary)  CL  from  Eq.  (5.1)  becomes  equal 
to  the  ratio  of  the  number  of  equal  time  intervals;  namely 
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Figure  5.1.  Complexit)  as  a function  of  LASA  mb  for  all 
reported  shallow  earthquakes  (O's)and  pre- 
sumed explosions  (X's)  in  the  data  base. 
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6/1  ior  'he  0 to  5 sor  'i  i ul  5 to  36  si'c  i inu*  intervals  ton 
side  red  here. 

Complexity  values  computed  for  the  two  large  presumed 
explosions  at  Novaya  Zemlya  (CL  = 2.4,  mb  = 6.1;  CL  = 2.7, 
mb  = 6.4)  are  greater  than  values  for  explosions  of  comparabl 
magnitude  (Fig.  5.1)  detonated  in  different  geographical  re- 
gions. The  relatively  complex  character  of  P-wave  trains 
from  events  at  Novaya  Zemlya  has  been  studied  by  other  in- 
vestigators (Douglas,  et  al^,  1973;  Frasier  and  Yang,  1975). 
The  explanation  that  appears  most  consistent  with  all  the 
available  data  is  that  body  waves  are  scattered  in  the  source 
region  and  travel  relatively  high-Q  paths  to  receivers  as 
compared  to  the  direct  P-wave.  Whether  or  not  this  explana- 
tion is  correct,  the  fact  is  that  the  event  population  would 
have  tc  be  substantially  regionalized  for  complexity  versus 
mb  to  be  an  effective  discriminant. 

5.3  SPECTRAL  RATIO 

Kelly  (1968)  and  Lacoss  (1969)  evaluated  a spectral 
ratio  criteria  that  discriminates  between  earthquakes  and 
explosions  based  on  recordings  of  short  period  P-waves. 

The  definition  of  spectral  ratio  employed  by  these  investi- 
gators is. 


(5.2) 


where  A(f^)  are  the  spectral  estimates  at  frequencies  f^ 
based  on  10  to  20  seconds  of  data  beginning  with  the  onset 
of  the  P-wave.  In  practice  the  numerator  is  summed  over 
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quency  components.  Several  combinations  of  frequency  para 
meters  (?,  m;  p,  q)  and  two  time  windows  (12.8  and  U . 4 secs) 
were  tested  using  the  LASA  data  base.  A sampling  of  the  re- 
sults, including  those  combinations  that  yielded  the  best 
discrimination , are  shown  in  Figs.  5.2a-c.  The  frequency 
bands  for  the  results  in  Fig.  5.2a  are  similar  to  those  used 
by  Lacoss  (1969),  while  the  data  in  Figs.  5.2b  and  c shew 
the  results  of  emphasizing  the  high  frequency  (>  2 Hz)  por- 
tion of  the  event  spectra. 

Examination  of  Figs.  5a-c  indicates  that,  as  in  the 
case  of  the  VFM  discriminant,  as  the  separation  of  the  high 
and  low  frequency  bands  increases,  discrimination  is  enhanced. 
The  degiee  of  discrimination  achieved  in  the  case  of  Fig. 

5.2c  is  in  fact  comparable  to  that  achieved  with  the  VFM 
technique  (Fig.  4.2).  This  is  not  terribly  surprising  since 
the  VFM  discriminant  can  be  considered  as  a finely  tuned 
( f requcncy - wi se ) spectral  ratio  tec  hnique . The  major  ad van  - 
tage  of  the  VFM  over  the  spectral  ratio  discriminant  is  the 
fact  that  the  VFM  technique,  in  addition  to  providing  for 
event  identification,  can  be  used  for  determinat i on  of  explo- 
sion yields.  The  large  amount  of  scatter  of  the  explosion 
populations  in  the  case  of  spectral  ratio  versus  mb  (Figs. 

5- 2a-c)  with,  if  anything,  an  inverse  relationship  between 
spectral  ratio  and  event  size,  would  seem  to  preclude  the 
possibility  of  using  this  technique  as  an  explosion  yield 
indicator . 

5.4  HIGHER  MOMENTS  OF  FRFQUF.NCY 

Recently  Canadian  researchers  (Anglin,  1971;  Weiehrrt  , 
1971;  Manchee,  1972)  proposed  combinations  of  time  and  fre 
quency  domain  algorithms  as  short  period  discriminants  between 
earthquakes  and  underground  explosions.  These  algorighms  are 
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Figm-e  5.2a.  Spectral  ratio  versus  LASA  mb  for  the  Eurasian 
? ,Ja,ibaSe;  Shallow  earthquakes  are  indicated  by  O's  and  pre- 

25  fO  98PHj?10SS-bs*Xri*A«:I3  ?his.case  1 = 12  f°-47  Hz) , m = 
10.38  Hz),  p * 3/  Cl. 45  Hz)  and  q = 50  (1.95  Hz). 
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rigure  5.2b.  Spectral  ratio  versus  LASA  mu . In  this  test 
1 = 8 CO . 31  Hz),  m = 19  (0.74  Hz),  p - 37  (1.45  Hz)  and 
q - 64  (2.50  Hz). 


115 


SPECTRAL  PATIO  (10 


t 


R- 2646 


?.o* 


1.5 


1.C+ 


LASA  DATA  WINDOW  ■ 12.80  SECONDS  PLOT 


x x * 
x 


.51 


.ol 

3.5 


X 

0 

0 

0 0 


0 

0 

0 8 


•*.5 


0 

8 o 

0 

0 


0 

0 0 
0 


0 o 


5.5 


0 

0 


• :« 
0 


MB  (10 


o) 


Figure  S . 2c . Spectral  ratio  versus  LASA  mh.  In  this  test 
q ; ?6C?2397HHz5.m  = 19  (°'74  H,)  ’ P = 51  ^-99  Hz)  and 
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referred 
modified 
given  by 


lO  as  the  third  moment  of  frequency  (TMF)  and 
moment  of  frequency  (MMF)  . The  definition  of  MM1: 
Manchee  includes  TMF  and  is 


MMF 


^ C[A(fk)  - N C f k ) ] * f 


k=r 


X tA(v  • N(fk’] 

k = r 


where  A(f^)  and  N(f^)  are  the  amplitudes  of  the  kth  spec- 
tral estimates  of  the  signal  and  preceding  noise  sample, 
respectively;  J,  M are  the  orders  of  the  moments  being  con- 
sidered in  their  respective  terms;  r,  p correspond  to  the 
lowest  frequency  components;  and  s,  q to  the  highest  fre- 
quency components.  When  used  as  a discriminant  MMF,  or  TMF 
given  by  the  first  term  in  Eq.  (5.3),  is  plotted  against 
either  the  body  wave  magnitude  or  complexity  of  an  event. 

The  rationale  behind  these  discriminants  is  based  on 
the  recognized  spectral  differences  between  earthquakes  and 

explosions.  The  effects  of  the  multipliers  f^  and 

M ^ 

(fq  - f^  + f ) in  Eq.  (5.3)  are  to  enhance  the  spectral 

differences  at  high  and  low  frequencies,  respectively,  be- 
tween earthquakes  and  explosions  while  minimizing  the  central 
portion  of  the  event  spectra. 

The  TMF  and  MMF  algorithms  were  coded  and  applied 
to  the  LASA  data  base.  Two  time  windows,  6.4  seconds  and 
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and  12.8  seconds,  and  various  combinations  of  frequency  para 
meters  were  tested.  As  explained  previously  no  noise  cor- 
rections were  made  in  these  tests.  The  results  are  shown  in 
Figs.  5.3a  and  b and  Figs.  5.4a  and  b. 

TMF  versus  m b and  linear  complexity,  C^,  are  plotted 
in  Figs.  5.3a  and  b,  respectively,  for  the  Furasian  popula- 
tion of  shallow  earthquakes  (0's)  and  presumed  explosions 
(X's).  The  results  of  these  two  figures  represent  the  best 
discrimination  achieved  for  many  different  sets  of  para- 
meters (r,  s;  Eq.  (5.3))  tested. 

In  the  case  o^  TMF  versus  m^  (Fig.  5.3a)  the  event 

populations  exhibit  complete  separation  of  m,  ' s > 4.7.  The 

b 

small  magnitude  explosions  that  plot  within  the  earthquake 
population  also  failed,  because  of  background  noise,  to  dis- 
criminate using  the  VFM  technique  (Fig.  4.2).  For  TMF  versus 
CFig.  5.3b)  the  two  event  populations  nearly  separate  over 
the  entire  magnitude  range.  The  degree  of  separation  using 
the  TMF  technique,  however,  is  not  as  great  in  either  case 
(Figs.  5.3a  and  b)  as  the  separation  of  the  event  populations 
realized  with  the  VFM  discriminant. 

Values  of  MMF  were  computed  for  the  Furasian  event 
population  using  Eq.  (5.5)  and  are  plotted  in  Figs.  5.4a 
and  b versus  and  linear  complexity,  (Cj ) respectively. 
While  there  is  substantial  overlap  of  the  event  populations 
in  Fig.  5.4b  (MMF  versus  C^) , the  data  in  Fig.  5.4a  indicate 
that  MMF  versus  m^  is  a very  effective  discriminant.  In 
fact  the  degree  of  separation  of  the  earthquake  and  explosion 
populations  in  Fig.  5.4a  is  as  good  as  or  better  than  that 
achieved  with  the  VFM  technique. 

Both  Oi.  the  higher  moments  of  frequency  techniques 
CTMF  and  MMF) , while  useful  as  positive  discriminants,  do 
not  appear  to  be  dependent  on  event  size,  except  in  an  in- 
verse manner.  Thus,  the  combined  application  of  the  higher 
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*T™  h3a'  BIF  versus  for  the  Eurasian  population  of  shallow  earthquakes 
CO  s)  and  presumed  explosions  (X's).  The  parameters  used  in  this  test  are 
r = 8 (.31  Hz)  and  s = 64  (2.5  Hz)  (see  Fq.  (5.3)). 


UJA  C*l«  WINDOW  . It. (•  SECONDS 


CINt**  COBPlCXITT  (It  ** 


Figur-1  5.5b.  IMF  versus  linear  complexity.  Parameters  used  are  r = 12 
(.47  Hz)  and  s = 64  (2.5  Hz). 
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lAJA  DATA  WINDOW  ■ I?. BO  SECONDS 


5*4a*  ^ versus  % for  the  Eurasian  population  of  earthquakes 
(0  s)  and  explosions  (X's).  The  parameters  used  in  this  test  are  r = 
12  (.47  Hz),  s = 64  (2.5  Hz),  p - 12  (.47  Hz)  and  q » 26  (1.02  Hz). 
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figure  5.4b.  MMF  versus  linear  complexity  for  Eurasian  event  population. 


120 


R-2646^ 


moment  discriminants  and  the  VFM  discriminant  would  offer  n 
very  powerful  approach  to  the  problems  of  event  ident  ificat  ion 
and  explosion  yield  determination. 

5.5  M - m, 
s b 

Lacoss  (1969)  applied  the  Mg  - m^  discriminant  to  the 
Eurasian  population  of  events  under  consideration  here.  As 
in  the  case  for  all  the  short  period  techniques,  many  of  the 
deep  earthquakes  were  found  to  appear  explosion-like  by  the 

Ms  ’ mb  criterion-  Thus,  the  results  shown  in  Fig.  5.5  are 
for  shallow  events  only. 

While  the  majority  of  earthquakes  separate  comfortably 
from  the  explosion  population  in  Fig.  5.5,  one  fairly  large 
magnitude  earthquake,  a reported  shallow  event  in  the  Kurile 
Islands,  lies  within  the  explosion  population.  In  addition, 
a rather  large  magnitude  (m^  * 5.5)  presumed  explosion  plots 
quite  close  to  the  limit  of  the  earthquake  population.  More 
important  than  this  overlap  of  populations,  however,  is  the 
fact  that  the  Mg  - criterion  is  limited  in  its  application 
to  explosions  of  relatively  large  magnitude.  In  this  particu- 
lar case  (Fig.  5.5)  there  are  no  detected  surface  waves  from 
explosions  of  mfe  < 5.1  and  several  instances  of  no  observed 
surface  waves  from  explosions  with  m^'s  in  the  range  5.4  to 
5.8.  On  this  count  the  VFM  and  other  short  period  discriminants 
previously  described  appear  quite  superior  in  that  their  range 
of  applicability  extends  to  events  approximately  one  full 

magnitude  unit  smaller  than  in  the  case  of  M - m 

s b * 

5-6  MULTIPLE  EXPLOSION  SCENARIO 

ihe  various  short  period  discriminants  (complexity, 
spectral  ratio,  TMF  and  MMF)  were  applied  to  the  multiple 
explosion  scenario  previously  described  in  Section  4.7  of 

this  report.  The  results  of  this  test  are  given  in  Table 
5.1. 
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LA$A 

Figure  5.5.  LASA  Ms  versus  LASA  m^,  for  events  with  detected 
surface  waves  (events  with  USCGS  depth  _>  100  km 
excluded).  Arrows  pointing  down  attached  to  the 
horizontal  bars  denote  no  Rayleigh  waves  de- 
fected. 
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3.11 

0.84 
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3.27 

5 

11.97 

1.93 

5.74 

3.60 

The  parameters  used  in  this  test  were  the  same  as  those  used 
in  each  of  the  previous  tests.  For  instance,  the  TMF  dis- 
criminant was  applied  to  the  first  12.8  seconds  of  the  simu- 
lated multiple  time  series  with  r = 12  (0.47  Hz)  and  s = 64 
(2.5  Hz).  The  integers  (1-5)  in  the  far  left-hand  column 
refer  to  the  different  shot-array  to  receiver  azimuths. 

As  explained  in  Section  4.7,  the  scenario  was  designed 


to  produce  teleseismic  signatures  that  resemble  those  from 
earthquakes  when  examined  in  light  of  certain  proposed  dis- 
criminants. Comparing  the  multi-azimuth  values  of  linear 
complexity  foi  the  simulated  multiple  explosion  with  the 
value  for  the  original  presumed  explosion  in  Table  5.1,  we 
find  that  the  scenario  has,  at  least  in  the  case  of  this 
particular  discriminant,  produced  the  desired  result.  That 
is  the  multiple  time  series  appears  earthquake-like  and,  in 
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fact,  at  certain  azimuths  (3,5)  more  cart  hq  mko  like  than 
most  of  the  events  plotted  in  Fig.  5.1. 

The  situation  (Table  5.1)  in  the  case  of  the  other 
short  period  discriminants  is  quite  different,  however.  The 
values  of  spectral  ratio,  TMF  and  MMF  for  the  multiple  explo- 
sion at  all  five  azimuths  remain  within  the  respective  explo- 
sion populations  (Figs.  5.2a,  5.3a  and  5.4a).  Thus,  while 
this  scenario  may  appear  earthquake - 1 ike  in  terms  of  first 
motion,  complexity  and  Mg  - mb,  the  short-period  discriminants 
examined  here  would  correccly  identify  this  event  as  an  explo- 
sion. 


5.7  SUMMARY 

The  relative  effectiveness  of  the  VFM  and  several 
other  short-  and  long-period  discriminants  was  compared 
using  the  same  population  of  wel 1 -recorded  Eurasian  events. 
The  most  significant  results  obtained  are  the  following: 

1.  The  short -period  techniques  (VFM,  spectral 
ratio,  TMF  and  MMF)  can  be  applied  to  earth- 
quakes and  explosions  over  a much  greater 
magnitude  range  than  the  Ms  - mb  discriminant. 

In  particular,  the  short-period  discriminants 
provide  positive  event  identification  for  explo- 
sions approximately  1 magnitude  unit  smaller 
than  can  be  tested  with  the  Ms  - mb  technique. 

2.  The  combined  application  of  the  VFM  and  MMF 
discriminants  should  prove  to  be  a very  ef- 
fective method  for  event  discrimination  and 
explosion  yield  determination. 

3.  Short-period  discriminants  based  on  spectral 
differences  between  earthquakes  and  explosions 
can  be  used  to  discriminate  earthquakes  and 
multiple  explosions  as  well. 
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VI.  DATA  ACQUISITION  AND  IhlTlNO 

From  the  outset  of  this  research  program,  the  acquisi- 
tion of  data  has  been  guided  by  the  requirement  for  a compre- 
hensive data  base  for  several  criteria  including:  testing 

and  verification  o^  our  waveform  modeling  capabilities;  the 
applicability  of  new  and  existing  discriminants;  and  the 
specification  of  a matched  filter  for  improved  event  detec- 
tion and  identification.  To  date  we  have  acquired  both  digi- 
tal and  analog  data  for  more  than  400  events  originating  in 
Eurasia  and  North  America  and  recorded  at  several  different 
world -wide  locations.  The  present  in-house  library  of  events 
consists  of  the  following: 

1.  Best  beam  short-period  (P-waves)  digital 
recordings  from  LASA  of  190  Eurasian  events 
(155  earthquakes,  35  explosions). 

2.  Best  beam  short-period  (P-wave)  digital  re- 
cordings from  the  original  Oyer  subarray  in 
Norway  of  44  Eurasian  events  (36  earthquakes, 

8 explosions) . 

3.  Three  component  short-  and  long-period  digital 
recordings  from  numerous  LRSM  stations  for  21 
NTS  explosions,  the  Aleutian  explosion  Longshot, 
and  a presumed  Eurasian  explosion. 

4.  Copies  of  seismograms  from  the  Palmer  network 
of  short-period  stations  located  throughout 
Alaska  and  the  Aleutians  for  more  than  100 
Eurasian  and  North  American  earthquakes  and 
explosions . 

5.  Copies  of  seismograms  from  selected  WWNSS  sta- 
tions for  a suite  of  NTS  explosions. 

6.  Copies  of  high-gain  long-period  seismograms  of 
NTS  and  Eurasian  explosions. 
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Digital  data  from  the  Yellowknife  array  in 
Canada  for  a suite  of  NTS  explosions  and 
earthquakes  located  in  western  North  America. 

While  most  of  the  above  data  are  for  events  recorded 
at  teleseismic  distances,  we  expect  to  receive  digital  data 
(courtesy  of  Lane  Johnson)  for  NTS  explosions  and  nearby 
earthquakes  recorded  at  Jamestown,  California  (JAS) , as 
well  as  WWNSS  seismograms  for  stations  closer  than  about  20° 
to  certain  Eurasian  source  regions.  These  data  will  provide 

important  checks  on  our  modeling  capabilities  at  regional 
distances . 

An  editing  program  for  handling  variable  formatted 
digital  data  has  been  written.  This  program  converts  the 
tape  data  that  comes  into  S3  to  a format  that  is  compatible 
with  the  various  on-line  seismic  analysis  codes.  Plotting 
routines  are  included  in  the  editing  program  and  allow  us 
to  check  the  quality  of  the  digital  data  before  processing. 
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APPENDIX  A 

THE  REFLECTION  OF  SPHERICAL  WAVES  AT  A FREE  SURFACE 


It  is  oi  ten  useful  and  enlightening  to  represent  an 
explosion  or  earthquake  in  terms  of  an  equivalent  elastic 
point  source.  For  studies  of  teleseismic  body  and  surface 
waves,  this  representation  is  computationally  convenient  and 
nearly  always  a reasonable  approximation  of  reality. 

Among  the  most  important  phases  in  the  teleseismic 
body  wave  signature  of  a complex  event  are  the  free  surface 
reflected  phases  pP,  pS,  sS  and  sP.  Using  frequency  inde- 
pendent free  surface  reflection  coefficients,  these  phases 
may  be  related  to  the  direct  P and  S waves  emanating  from 
the  source  region,  thereby  deducing  a great  deal  of  informa- 
tion about  the  source. 

ihe  reflection  coefficients  for  plane  waves  impinging 
on  a free  surface  may  be  found  in  many  textbooks  (e.g., 

Ewing,  Jardetsky  and  Press  (1957),  Chapter  2).  However,  we 
have  not  seen  in  the  literature  a straightforward  presenta- 
tion o f ^ the  reflection  coefficients  for  the  case  of  a point 
source.  Further,  we  believe  there  may  be  a certain  amount 
of  confusion  on  this  point. 

In  sophisticated  methods  employing  generalized  rays 
(e.g.,  Wiggins  and  Helmbergcr  (1974)),  or  propagator  matrices 
(e.g.,  Fuchs  (1966)),  the  sphericity  of  the  wave  front  is 
automatically  included  in  the  computation.  However,  it  is 
oxten  desirable  to  treat  the  source  region  as  a homogeneous 
half  space  and  compute  the  reflected  phases  directly  from 
the  upgoing  P and  S waves  generated  by  the  source.  In  these 


It  is  true  that  these 
with  a fair  amount  of 
of  Gupta  (1967). 


coefficients  can  be  deduced,  albeit 
algebra,  from  results  such  as  those 
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cases  the  plane  wave  coei' f icicnts  suffice  for  the  pi’  and 
sS  phases,  hut  a sphericity  correction  is  needed  to  compute 
pS  and  si’.  In  the  following  paragraphs  the  necessary 
sphericity  corrections  for  obtaining  the  plane  wave  phases 
pS  and  si’  from  the  point  source  generated  spherical  waves 
incident  on  a free  surface  will  be  given.  While  there  are 
many  ways  to  derive  these  correction  factois,  a method 
employing  a saddle  point  approximation  to  an  appropriate 
branch  line  integral  is  convenient  and  illuminating. 

Plane  Wave  Reflection  Coefficients 

Consider  a homogeneous,  isotropic  half  space  with 
the  free  surface  at  z = 0 in  a cylindrical  coordinate  system 
(r, 0 , z) . Assuming  azimuthal  symmetry,  the  Fourier  trans- 
formed displacements  in  the  r and  z directions  may  be 
written  in  terms  of  potentials  as  follows: 


q(r>z,w)  ■ !£  * slfj  . 
w(r,z,w)  = + k2^  , 


where  k = w/c,  with  c the  usual  horizontal  phase  velocity. 

For  plane  waves  incident  at  a free  surface,  the  reflec- 
tion coefficients  are  easily  derived  using  these  potentials. 

In  terms  of  displacements,  these  are 


IV1  ‘ A|V  ■ 

'ussl=Alu5l  . 


(A. 2) 
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where 


A = ' (r8  ~ 

4VS  * (rB  ‘ D2 

and 
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ps 


r a 

"TP  Blupl  • 


r«s 

UsP  = -5-  Blusl  • 


where 


4(r2  - i) 
B = S 


4VB  + (rB  • 


(A. 3) 


In  these  equations  ra  = (c2/a2  - 1)1/2  and  rg  = (cz/B2  - 1)1/2 
while  i Up | , |Upp| , etc.,  are  the  amplitudes  of  the  incident 
and  reflected  plane  waves  at  the  free  surface. 


Sphericity  Corrections 

Consider  now  a point  source  at  a depth  z = h in  the 
half  space.  At  teleseismic  distances  the  waves  from  the 
source  and  nearby  free  surface  may  be  treated  as  plane  waves. 
However,  the  sphericity  of  the  wave  front  should  be  taken  into 
account  when  computing  the  waves  reflected  from  the  free  sur- 
face. For  the  reflected  phases  in  which  the  wave  type  remains 
unchanged  (Upp  and  U^g) , the  plane  wave  coefficients  of 

are  equally  valid  for  a point  source.  However,  for  the 
converted  phases  (Usp  and  UpS)  a sphericity  correction  is 
required.  .his  correction  takes  the  simple  form  (r  /r  )— i. 
ihat  js,  fur  spherical  waves  incident  on  a free  surface  (A. 5) 
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where  L(  is  a path  in  the  fourth  quadrant  of  the  complex 
C plane.  This  representation  is  for  body  waves  only  ami 
is  valid  in  the  far  field  (hr  >>  1). 

The  integral  (A. 6)  has  a saddle  point  and  may  be  ap- 
proximated using  the  method  of  steepest  descent.  The  perti- 
nent geometry  is  shown  in  rig.  A.l. 

The  single  real  saddle  point  of  (A. 6)  may  be  written 

* sint,p  • (A. 7) 

where  £ is  the  branch  point  associated  with  r . Taking 

vA  CX 

the  path  of  steepest  descent  through  the  saddle  point,  it  is 
found  that  (returning  to  real  wave  number) 


B ” 1 (kaR  + M ) 
V(r,z,w)  = j— e 1 8 2 

o 


(A. 8) 


where  k = m/a,  k = m/B  and  k = sin0  k = sin0  k„. 
a B o pa  s 3 

lengths  R , R are  defined  in  Fig.  A.l.  Further. 

12  ’ 


The 


D2  = 


COS20  COS20 

s p 


(A.  9) 


The  meaning  of  (A. 8)  is  that  for  a given  z,  h and  r 
there  is  only  one  ray,  specified  by  0 , which  travels  between 
source  and  receiver.  The  phase  given  in  (A. 8)  is  the  exact 
phase  travel  time  for  the  pS  ray  while  the  amplitude  is  ap- 
proximate. 
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Figure  A.l  Geometry  for  the  S waves  arriving  a 
depth  z and  a distance  r from  a 
source  of  P waves  at  a depth  h. 
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Plane  Wave  Interpretation 


l'he  plane  wave  approximation  is  valid  when  the  point 
ol  observation  is  removed  a great  distance  1 rom  the  source. 
Ior  (A. 8)  this  approximation  is  of  interest  for  two  limiting 


cases : 


Case  1 r >>  r 
1 2 


Since  (A. 8)  is  va  _d  only  for  kr  >>  1,  this  case  may 
be  viewed  as  equivalent  to  h >>  2;  that  is,  a deep  source  as 
viewed  near  the  surface.  Since  the  total  displacement  is 


upS  = ‘1  cos9s  " w i',in6s » it  can  be  shown  that 


— fll  k ~ i ( k R + k _ R ) 

Y - ‘Vc  {-  « “ ‘ 3 ! 

1 


(A. 10) 


where  the  superscript  indicates  that  the  displacement  applies 
for  Case  1. 


Cased  r >>r 
2 i 


In  this  case 


U(2)„.k  r — 
pS  Ys  R 


(A. 11) 


ihe  j.atter  case  is  equivalent  to  z >>  h and  is  appro- 
priate for  teleseismic  observations  of  nearly  any  earthquake 
or  explosion  source. 

From  (A. 10)  it  is  easily  shown  that  the  reflection  co- 
efficient at  the  free  surface  for  the  pS  phase  is  identical 
to  that  given  in  (A. 3),  that  is,  the  plane  wave  coefficient. 


Cn  the  other  hand,  Case  2 gives  the  coefficient  in- 
cluding the  correction  for  the  sphericity  of  the  wave  front 
incident  on  the  free  surface.  For  the  pS  phase  this  correc- 


tion is  r_5/ra. 
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Ci^ri  > ing  out  a similar  analysis  for  a center  of  rota- 
tion, it  may  be  shown  that  the  sphericity  correction  to  be 
applied  to  the  plane  wave  coefficients  for  the  sP  phase  is 
W Thus’  the  freu  surface  reflection  coefficients  for 
the  pS  and  sP  phases,  including  the  effect  of  sphericity  of 
the  wave  front,  are  given  by  (A. 4). 

Concluding  Remarks 

It  is  commonly  observed  that  plane  wave  approximations 
are  valid  for  sources  viewed  at  a great  distance  and,  there- 
fore, i or  teleseismic  observations  of  almost  any  earthquake 
or  explosion  source.  Elementary  calculations  of  the  free 
surface  reflected  phases  made  by  applying  a simple  multiplica- 
Live  constant  to  the  direct  source  generated  wave  are  often 
useful.  For  the  phases  pP  and  sS  the  plane  wave  reflection 
coefficients  will  suffice.  However,  for  the  converted  phases 
PS  and  sP,  a sphericity  correction  (yyl1  must  be  included. 
Failure  to  include  this  factor  will  result  in  an  error  which 
is  on  the  order  of  (2)-1. 
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APPENDIX  B 

THE  SEISMIC  BODY  WAVES  GENERATED  BY  A 
POINT  SOURCE  IN  A LAYERED  MEDIUM 


I.  INTRODUCTION 


Details  of  the  crustal  structure  in  the  vicinity  of 
both  the  source  and  the  receiver  can  have  an  important  effect 
on  the  recorded  body  waves.  In  the  receiver  region  one  can 
use  the  theory  of  Haskell  (1962)  in  which  the  emergent  waves 
trom  the  upper  mantle  arc  treated  as  plane  waves  incident 
from  below  on  a stack  of  plane  layers  representing  the  crust. 
This  procedure  has  been  in  widespread  use  for  a number  of 
years  and  is  quite  successful  in  regions  where  the  plane, 
parallel  layer  assumption  is  not  too  far  removed  from 
reality. 

Assuming  a source  region  crustal  model  made  up  of  a 
stack  of  plane  elastic  layers,  a theory  analogous  to  that  for 
the  receiver  crust  can  be  constructed.  Due  to  the  presence 
of  the  source  in  one  of  the  layers,  this  theory  is  a good 
deal  more  complex.  Fuchs  (1966)  developed  a transfer  function 
representing  the  effect  of  the  layered  crustal  model  on  the 
P waves  at  large  distances  for  three  types  of  point  sources: 
a center  of  dilatation,  a single  couple  and  a double  couple 
of  arbitrary  orientation.  Kogeus  (1968)  used  Fuchs'  theory 
to  compute  teleseismic  synthetic  seismograms  for  explosion 
sources.  Hudson  (1969a, b)  extended  Fuchs'  results  to  apply 
to  quite  general  sources  of  finite  extent.  However,  the 
representation  of  complicated  sources  at  arbitrary  orientation 
in  a form  compatible  with  the  layered  crust  theory  given  by 
Hudson  appears  to  be  a complex  and  difficult  task.  Thus,  the 
usefulness  of  Hudson's  theory,  especially  for  routine  numeri- 
cal computations,  seems  to  be  limited  by  the  inherent  com- 
plexity of  the  source  representation. 
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Archambcau  (1964,  1968)  has  pointed  out  that  the 
radiatio.  rield  external  to  any  kind  of  volume  source  can  be 
represented  in  an  expansion  in  spherical  harmonics.  The  ex- 
pansion coefficients,  called  multipole  coefficients,  then 
provide  a unique  and  compact  numerical  representation  of  a 
source  of  arbitrary  complexity  and  orientation.  This  repre- 
sentation, called  an  equivalent  elastic  source  representation, 
provides  a convenient  link  with  powerful  elastic  wave  propa- 
gation techniques.  The  way  in  which  the  equivalent  elastic 
source  may  be  obtained  for  any  source  for  which  the  outgoing 
elastic  wave  can  be  computed  is  discussed  in  Section  II. 

Harkrider  and  Archambeau  (1975)  computed  the  surface 
waves  excited  by  an  equivalent  elastic  source  in  the  form  of 
Section  II  embedded  in  a stack  of  plane  layers.  The  corre- 
sponding theory  for  far  field  body  waves  is  presented  here. 
Having  the  body  waves  exciting  the  plane  layered  crustal  model, 
one  can  use  generalized  rays  (e.g.,  Wiggins  and  Helmberger 
(1974))  or  other  means  to  propagate  the  waves  through  the 
upper  mantle.  Then  one  may  apply  appropriate  transfer 
functions  for  the  receiver  crustal  model  and  sensing  instru- 
ment to  compute  synthetic  seismograms.  In  short,  this  theory 
provides  a merging  of  arbitrary  nonlinear  source  calculations 
with  the  linear  elastic  wave  propagation  techniques  developed 
in  seismology. 

It  should  be  pointed  out  that  the  results  to  be  pre- 
sented here  are  completely  equivalent  to  those  of  Hudson. 

The  difference  is,  of  course,  the  specification  of  the  source 
in  the  computationally  convenient  multipolar  coefficient  form. 


I 
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x1,  EQUIVALENT  elastic  source  representation 

The  Fourier  transformed  equations  of  motion  in  a homo- 
geneous, isotropic,  linearly  elastic  medium  may  be  written 


— I Vy(4)  + / L_ 


V X x 


(2.1) 


-here  u is  particle  displacement  and  ka  and  ka  are  the 

congressional  and  shear  wave  numbers.  The  Cartesian  poten- 
tials  J ' ~ 


% are  defined  by 


X(4)  = V • u 


X = J v X u 


(2.2) 


and  may  be  easily  shown  to  satisfy  the  wave  equation 

* k|7CJ)  - 0 . j . 1,  2,  3,  4,  (2.3) 

where  k = k^  = u/a  , q s kg  . u/e  , i - 1,2,3.  This  equa- 
tion has  as  a solution  the  following  expansion  in  spherical 
eigentunctions  (e.g.,  Morse  and  Feshbach  (1953)), 


- £><2>w  ±1 41 

£-0  m=o  ' 

+ Blm^  sin  Pj(cos6) 


(u)  cos  m<f> 


(2.4) 


whore  the  h£  ^ are  spherical  Hankel  functions  of  the  second 
Kind  and  the  P&  are  associated  Legendre  functions.  The 
vector  R has  as  components  the  spherical  coordinates  R,e,<p. 

Equations  (2.4)  together  with  (2.1)  provide  an  equiva- 
lent elastic  source  representation  of  the  outgoing  elastic 
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Wiivcs.  ihe  values  of  the  multipole  coefficients 


B 


Jim 


( ‘ » J 1 >2,3,4,  prescribe  the  displacement  field  at  all 


points  in  the  homogeneous  elastic  medium  where  (2.1)  applies. 

The  point  source  representation,  liq.  (2.4),  provides  a 
generalization  of  the  commonly  used  center  of  dilatation 
(£  = 0),  couple  (Jl  = 1)  and  double  couple  (i  - 2)  point 
sources,  as  well  as  more  complex  sources.  lor  example,  for 
a horizontal  double  couple  source  of  unit  amplitude,  the  non- 
zero coefficients  in  (2.4)  are 


^ - (SVa’) 


2 1 


2 2 


2 2 


Given  the  displacement,  u^,  or,  alternatively,  the 
potentials  , one  may  determine  the  multipole  coefficients 


using  (2.4).  For  example,  using  the  orthogonality  of  the 
spherical  eigenfunctions,  one  may  derive 


(j) 


Urn 


(j 


Jim 


w)  j 


2tt 


''Ifcosa, 


(cos  mb) 

J sinOdedb 
(sin  mb| 


(2.5) 


where 


r a (21+1)  (Jl-m) ! 
£m  27T(?.+m) 


m f 0 


C = (2Jl+l)/47r 


That  is  > ^given  the  displacements,  u(R,u)  on  a sphere  of 
radius  R,  one  can  use  (2.2),  (2.5)  and  (2.4)  to  compute  the 
displacements  at  any  point  in  the  elastic  regime.  This  pro- 
cedure was  first  suggested  for  geophysical  applications  by 
Archambeau  (1968). 


Using  the  equivalent  elastic  source  representation, 
highly  nonlinear  finite  difference  calculations  of  the  near 
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source  earthquake  or  explosion  displacement  fields  may  be 
coupled  to  the  elastic  field  where  efficient  elastic  wave 
propagation  techniques  can  be  applied.  Carrying  the  non- 
linear computation  out  to  a radius  at  *which  the  material 
response  is  in  the  elastic  regime,  Eq.  (2.5)  may  be  applied 
to  compute  the  multipole  coefficients  and,  thereby,  the 
equivalent  elastic  source. 

For  a spherically  symmetric  explosion,  the  explosion 
induced  elastic  field  is  often  represented  by  the  reduced 
displacement  potential,  v(t),  defined  by 


Ur(r>t,  .iiii.i-fiiii  , 


where  t = t - R/a,  is  the  retarded  time.  Note  that  ¥(t) 
is  independent  of  R. 

Applying  the  Fourier  transform,  (2.6)  becomes 


_ (1  + ik  r)  -ik  R 

ur(R,u>)  = 2L_  e a 

r R2 


a . 


Then,  carrying  out  the  algebra  (observing  that  the  u of  (2.2) 
is  the  Cartesian  displacement  vector),  it  may  be  shown  that 
the  multipole  coefficients  reduce  to  a single  nonzero  term, 
the  monopole: 


\o  (“)  = * (2, 

For  more  complex  sources,  additional  nonzero  coeffi- 
cients are  present.  For  example,  in  the  case  of  an  explosion 
in^an  axisymmetric  (about  the  z axis)  cavity  the  coefficients 
A£0  » Bjj,1  > A1i>1  > 1 = 2,4,5,  ...  , are  also  nonzero  (Cherry, 

— ^ a -*• • > n)*  Also,  Archambeau  (1972)  has  given  a theory 

for  teutonic  stress  release  due  to  underground  explosions 

which  leads  to  nonzero  quadrupole  coefficients  A1^^ 

„ - - 2m  * 2m  ' 

m = 0,1,2. 


144 


I r. 


o 


K - 2 6 4 1> 


1 o thls  point  the  discussion  has  been  independent  of 
the  orienLation  of  the  source  with  reference  to  a free  sur- 
race, which  must  be  included  to  study  the  source  as  a generator 
of  body  and  surface  waves.  Minster  (1974)  has  derived  a 
general  transformation  by  which  the  multipole  coefficients  in 
a standard  coordinate  system  (e.g.,  fixed  with  respect  to  the 
free  surface)  may  be  obtained  from  any  rotated  coordinate  sys- 
tem. Therefore,  the  multipole  coefficients  may  be  computed 
in  a convenient  source  related  system  and  then  rotated  to  a 
fixed  geographical  system. 


J 


K-2046 


in.  AN  LLASTODYNAMIC  S0URC1:  IN  A PLANh  LAYLRLI)  HALF  SPACh 

Ilarknder  and  Archamocau  (1975)  have  develo]>ed  the 
theory  for  computing  surface  waves  for  an  equivalent  elastic 
point  source,  given  by  Lq.  (2.4),  embedded  in  a sta^k  of  homo- 
geneous, isotropic  elastic  plane  layers.  This  required  the 
formulation  of  the  displacements  in  terms  of  certain  integrals 
over  wave  number.  The  surface  waves  are  then  given  by  the 
residue  contribution  to  these  integrals.  The  body  waves 
exiting  the  plane  layered  model  are  given  by  the  branch  line 
contributions  to  the  same  integrals.  The  essential  steps 
for  deriving  the  k integrals  in  a form  suitable  for  compu- 
tation of  body  waves  are  outlined  below.  For  the  most  part 
the  notation  is  that  of  Harkrider  (1964)  and  ilarkrider  and 
Archambeau  (1975) . 

Consider  a semi-infinite  elastic  medium  made  up  of  n 
parallel,  homogeneous,  isotropic  elastic  layers.  Number  the 
layers  from  1 at  the  free  surface  to  n for  the  underlying 
half  space.  Place  the  origin  of  a cylindrical  coordinate 
system  (r , <J> , z ) at  the  free  surface  and  denote  the  layer  inter- 
faces by  z.,  i = 1,2,  ...,  n-1.  This  geometry  is  depicted  in 
Fig.  i.  Let  (q^,  v ^ , w^)  be  the  components  of  the  Fourier 
transformed  displacements  in  the  (r,$,z)  directions  in  the 
ith  layer.  Then,  following  ilarkrider  (1964),  the  potentials 
^i  * ^ are  defined  by 


34>.  aV 
Mitr.P,:)  - J--  * 


3‘pj  x sa. 

3r3z  + F 9$ 


d2^.  dn. 


_ 1 u-r.  1 y y.  OH. 

vi(r-*-z)  915?  ' sF  • 


wi(F,c?,z)  = 


30-  S2v 


- + k*  i = 1,  2,  ...,  n, 


Receiver^ 


i 


♦ 

z 


i 


The  geometry  and  coordinate  system  for 
source  at  a depth  h in  a stack  of  n-1 
layers  underlain  by  a half  space. 


a point 
plane 
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We  will  subsequently  be  interested  in  an  equivalent 
elastic  point  source  at  a depth  z = h.  Let  the  source  layer 
be  denoted  by  a subscript  s;  that  is,  z.  < h < z . It  is 
necessary  to  express  the  cylindrical  potentials  , Ti  ) 

in  the  source  layer  in  terms  of  the  spherical  potentials 
(Xg  ) from  Eq . (2.4).  The  equivalence  is  given  by 
Harkrider  and  Archambeau  (1975)  and  is  as  follows: 


4> 


s E f l^m  cos  + \ sin  h 
m=0  */0  1 ) <* 


-ikra|z-h| 


Jm(Xr)dk  , 


Z 


*s  ■ E / jEm 

m = n n ' 


m=0  0 
l 

■ e/ 

m=0  0 ' 


cos  m<j>  + F sin  m<fc 
m r 


-ikr  | z-h | 


e 


Jm(Xr)dk  , 


C^3^cos  m<p  + D sin  md> , 
m m l r 


+ bi^sin  mt)!  — 


-ikr  | z-h | T ,,  . 
6 Jm(kr) 


(3.2) 


dk  , 


where 


L 


E AC4)P?  (£»)  , 

a H=m  v a ' 


G 


*.-7S  (0»-"[Sgn(h-Z„-«  Bf«P-  Ip.)  , 
a 1= m ' a ' 


krJ 

8 


(3.3) 


^ " jjj  E Ci)m-ntsgn(h-Z)]m^A^Pj(^)  , 


G 


crCj) 

m 


ft  E 0)m"n[sgn(h-z)]m+t  bOJp;  (£i)  , 

k8  £=m  VS' 


r 
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where  j - 1,2,5  in  the  definitions  of'  j)(j) 

m ’ m ’ 


and 


k = w/c,  k^  = w/y  , 


1, 

'2 


ry  = (c*/Y2-  1)’  , y = a or  3 , 


c - horizontal  phase  velocity. 


(3.4) 


ihe  coefficients  E , Fm  in  the  expression  for  ip 
defined  from  the  Dp  ^ as  follows: 


<Ps  are 


, F = 0 


2kE  = 

o 1 

2kE  = - D(1^  - 2C(2) 

12  2 0 

2kF  = C(1)  + D^2-^  + 

i 

0 for  m > i, , 


(3.5) 


12  2 0 
and,  taking  = D ^ = 


m 


m 


2kE  = C(2J  - - D O) 

m m+1  m-1  m+1  m-1  » 

2kF  = C(i)  + H1)  + n(23  . p(2)  , „ 

m m+1  m-1  m + 1 um-l  1 2 .1  m ^ • 

Since  the  location  of  material  boundaries  depends  only 
on  z,  the  dependence  on  r and  4>  will  be  everywhere  the 
same  as  in  the  source  layer.  Therefore,  separate  the  poten- 
tials in  the  layers  as  follows: 


Vr’«’z))  i 

00  1 

'*P°  (*,z) 

^(r.b.z)  J = 

/ j 

. (b,z) 

(r  > 4>  > z)  ) m=0 

0 i 

^i(m)(4,z) 

2^  • • • y n • 
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The  potentials  satisfy  wave  equa- 

tions fo.  which  the  general  solutions  may  be  written: 

/-m-*  -ikr  z ikr  z 

#i  ^(4>,z)  = A-  e ai  + e ai 


♦ ^(♦,2)  - e’lkrBiZ  ♦ Sf)'  elk"BiZ 


(3.7) 


Si}"1'1  (4>.z)  = £• 


-(m)  -ikre, 2 -Cm)'  ikrB-z 

= f>''e  l + e > ' e l 


in  Eq,  (3.7)  and  subsequent  equations  subscripts  i,  i 
n,  denote  quantities  in  the  ith  layer. 

Then,  define 


= 1,2, 


(m) 


‘•(Sr)*  •‘1‘Vi’1  - 


S?»)  = At!  e'ikr6iZi-l  ~Cm) 

1 Yi 


?<■)  = k e'^V1-1  ;w  . 


(3.8) 


We  will  subsequently  be  interested  in  the  values  of  the  poten- 
tial in  the  half  space  (i  = n) . Applying  the  radiation  condi- 
tion at  infinity,  (3.7)  and  (3.8)  give 
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4>Cm)(4>,  z)  = 


, -ikr  z \ 

i-  e an  4(m) 

k*  " 

a 

n 


iYn  /ikr6.2  .-0") 


n to. 


n 


(3.9) 


^m‘)  \ e Bn  , 


n 


where 


z = z - z 


n- 1 


y.  = 262/c2  . 
1 i 


Now,  combining  Eqs.  (3.6)  and  (3.9),  we  have 


(3.10) 


$n(r,<fr,z)  = 


t /(-  ■< 

m=0  0 ' Ka  ' 

n 


-ikr  z 

Cm)  e an  Jm(kr)dk  , 


^n(r  ,<f> , z)  = 


x,  00 

Zf 

m=0  *^0 


i 12 


a/'  x — ik  ~ z 

uAmj  e pn  J (kr)dk, 

2 i n m v J * 

3 ' (3.11) 


^n(r  »4> » z)  » 


X,  00 

Zf 

m=0  i) 


i -ikrQ  z 

F Sn  e " Jm0<r)dk 


m 


In  Eqs.  (3.11)  we  now  have  the  cylindrical  potentials 
*n»  ^n>  in  the  nth  layer  (half  space)  in  terms  of  a sum 

of  integrals  of  the  Sommerfeld  type.  The  coefficients 
An  * n * en  depend  on  azimuth,  <J> , as  well  as  wave 
number,  k,  and  include  the  modification  of  the  source 
generated  pulse  by  the  material  discontinuities  in  the 
layered  half  space.  Following  Harkrider  (1964,  Eqs.  (62) 
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and  (122)),  these  coefficients  are  solutions  of  the  matrix 
equation- 


2 Cm)' 
n 

/ 

uR  (0)/c‘ 
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lu  1 

m | 

jCm) 
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e v * 
n 

= JL 

vL  (0)/c 
i 

+ a:1 

6V 

m 
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/ 

0 

LS1 
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' 
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_ 
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where  the  J and  JL  are  given  by  Eqs.  (61)  and  (124)  of 
Harkrider  (1964).  The  matrix  AR  is  defined  by  Harkrider 
and  is  the  layer  product  matrix  wHch  gives  the  displacement- 
stress  vector  for  P - SV  motion  at  the  source  depth  in  terms 
of  the  displacement-stress  vector  at  the  surface.  Similarly, 

ALS1  is  the  transfer  matrix  for  the  displacement-stress 
vector  associated  with  SH  motion. 

The  source  terms  in  (3.12)  are  given  by  Harkrider  and 
Archambeau  (1975)  as  follows: 
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where 


6U 


m 


, U5 

6 1 — / cos  m<|>  + 6 
' ' m 


6W  = 
m 


./  ws 

6 I — I cos  m<|>  + 6 
\ ' m 


(S') 


sin  m|  , 


sin  mi})  , 


= ^°m  cos  m<^  + sxn  m<^  > 


6Xm  = 6xRm  cos  m<>  + 6TRm  sin  m<f>  , 


(3.13) 


6Vm  = Mr  cos  m<*>  + 6 

' ' m 


(H 


sin  , 
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6Ym  = 6TLm  cos  m<j>  + 6TLm  sin  m4>  » 
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2pc2k2  -yF  - ik(y-i)  — 


6TRm  “ 2pc2kZ  [-YB;  - ik(Y-l)  ~ 


Mc 


'*) 


£•(3)0 

i2k2  J2 , 

6 rB 

ff(3)o 

i2k2  JS , 

8 re 


4tL  * ' C<S>#  , 


<s t?  = - i2k2u  . r,  ,, 

Lm  g m (3.14 

The  quantities  r^,  rg,  p,  y , y,  kg  refer  to  the  layer  in  which 
the  source  occurs, with  p and  y being  the  density  and  shear 
modulus.  The  coefficients  3^,  Bm,  etc.  are  given  by  (3.3) 

with  the  following  modification.  The  series  are  separated  into 
two  parts; 

A = A6  + J° 

TTl  lO 


with  the  e superscript  denoting  a series  made  up  of  terms 
with  m + l even  and  the  o superscript  denoting  a similar 
series  from  terms  with  m + i odd. 

Equations  (3.12)  may  be  viewed  as  a set  of  simultaneous 
linear  algebraic  equations  and  solved  for  A^m^,  , £(m) 

in  terms  of  known  quantities.  The  result  is: 
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Equations  (3.16)  gives  us  ^m)  , S « , £(“)  in  terms 
Ot  the  known  multipolar  coefficients  defining  the  source 
and  the  layer  matrices  of  the  layered  half  space.  Then  (3.11), 
together  with  (3.1), give  closed  form  solutions  for  the  dis- 
placements in  the  nth  layer  or  half  space.  It  remains  to 
evaluate  the  integrals  in  (3.11)  to  extract  the  solution  for 
the  steeply  emergent  body  waves  which  are  of  interest  here. 

This  evaluation  is  the  subject  of  the  following  section. 
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IV.  COMPUTATION  OF  BODY  WAVES 


To  compute  the  displacement  potentials  in  the  half 
space,  (3.11),  it  is  necessary  to  evaluate  integrals  of  the 
form 


-ikr  z 

fm(k,o))  e Y Jm(kr)  dk  . (4.1) 

Y - % or  0n  , m » 0,  1,  2 

For  m = 0 the  dilatational  potential,  $ , was 
evaluated  at  large  distances  from  the  source  by  Fuchs  (1966) 
using  saddle  point  methods.  Hudson  (1969a, b)  encountered 
integrals  very  similar  to  those  in  (3.11)  when  solving  for 
the  body  waves  due  to  a point  source  in  a layered  media.  As 
pointed  out  in  the  Introduction,  Hudson's  solution  is  com- 
pletely equivalent  to  that  obtained  here,  with  the  difference 
being  in  the  specification  of  the  source. 

Hudson  solved  for  the  far  field  body  wave  contribution 
to  integrals  of  the  form  of  (4.1)  using  contour  integration 
in  the  complex  plane  and  approximating  the  branch  line  con- 
tribution using  the  saddle  point  method.  The  details  of  this 
integration  may  be  found  in  the  works  of  Fuchs  and  Hudson 
and  will  not  be  reproduced  here.  It  will  suffice  to  give  the 
results  which  are 


i = imtl 
m 


f(k,w)  r 


e 

R 


■ik  R 
Y 


(4.2) 


where  R‘  » r‘  + z2.  The  geometry  is  shown  in  Fig.  1.  The 
saddle  point  approximation  is  valid  as  long  as  R >>  r;  that 
is,  as  long  as  the  receiver  is  sufficiently  far  from  the 
base  of  the  stack  of  plane  layers.  Alternatively,  since 
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r/R  = sin6n,  this  approximation  is  valid  for  steeply  incident 
body  waves. 

With  the  solution  (4.2),  we  may  now  write  the  far 
field  body  wave  contribution  to  the  potentials  (3.11): 


*n(r.$.z) 


m=0  ' Ka  ' 
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i”*1  r 4 - - _ 
“n  n * 


-ik  R 
a (m)  e an  , 


■ ik„  R 


^nO,4>,z) 


m=0  \k6  / n 
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(4.3) 


^n(r,4>,z)  = 


im+l  r ^ R 

V i, r e 

m=0  n 


n 


The  displacements  in  the  half  space  are  given  in  the 
cylindrical  coordinate  system  by  applying  (3.1)  to  (4.3).  It 
is  more  convenient  to  deal  with  the  displacement  field  polar- 
ized into  its  P,  SV  and  SH  components.  In  the  far 
field,  these  are  given  by: 


_ 9$ 

Up  - sin6  *-2  + cose  v-2  , 

v n 3r  n 9z  * 
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(4.4) 
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1.  Compute  the  source  layer  potentials  $ , 
*s.  «s  from  (3. 2-3.5) . 

2.  Compute  the  displacement-stress  discon- 
tinuity quantities  6Um,  6Wn,  etc.,  from 
(3.13-3.14)  for  each  m. 


3.  Compute  the  layer  product  matrices  An 
R L Kq-. 

ALci  » J J • The  latter  two  matrices1 are 

independent  of  the  source,  while  the  first 

two  depend  only  on  source  depth. 


4. 


5. 


6. 


Compute  the  source  terms  Y , Z , X from 

m m m 

(3.18) . 

Compute  A^m),  <7‘,  e^m)  from  (3.16). 

Compute  the  displacements  at  a selected  R 
from  (4.5). 


t 
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APPENDIX  C 


A COMPARISON  OF  THE  SEISMIC  BODY  WAVES  GENERATED 
A POINT  COMPRESSIONAL  WAVE  SOURCE  AND  A 
CONCENTRATED  FORCE  AT  THE  FREE  SURFACE 


While  most  investigators  studying  the  teleseismic 
signature  of  underground  nuclear  explosions  have  been  pri- 
marily concerned  with  the  most  common  example,  fully  con- 
tained blasts,  another  class  of  explosion,  the  cratering 
sho*.  is  becoming  of  increasing  interest.  Before  embarking 
on  more  complex  and  sophisticated  investigations  of  the 
properties  of  the  teleseismic  signal  generated  by  cratering 
shots,  useful  guidance  can  be  obtained  by  comparing  crater- 
ing and  contained  explosions  using  very  simple  models.  In 
particular,  elastic  half  space  solutions  are  often  quite 
helpful  m understanding  the  gross  features  of  the  elastic 
waves  generated  by  seismic  sources. 

Two  elementary  point  source  representations;  (l)  A 
point  compressional  wave  source  at  shallow  depth,  and  (2) 

A concentrated  force  on  the  free  surface,  will  be  studied 
and  compared.  The  first  is  useful  for  representing  a con- 
fined explosion  while  the  second  provides  a simple  represen- 
tation of  an  explosion  on  the  surface.  The  latter,  in  turn, 
is  expected  to  give  some  indication  of  the  features  which 
make  cratering  and  contained  explosions  different. 

The  elastic  waves  in  a half  space  from  the  two 
sources  will  be  given  in  elementary  form  and  compared. 

Most  of  these  results  follow  directly  from  material  pre- 
sented in  Ewing,  Jardetsky  and  Press  (1957),  Chapter  2. 

The  important  results  will  then  be  cast  in  terms  of  the 
familiar  reduced  displacement  potential  to  facilitate  com- 
parison. Note  that  aU  secular  quantities  are  given  in  the 
Fourier  transform  domain.  For  convenience,  the  bar  gen- 
erally used  to  indicate  Fourier  transformation  is  omitted. 
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Consider  a cylindrical  coordinate  system  (r.0,z). 
For  the  half  space  problems  let  the  origin  be  at  the  sur- 
face with  z positive  downward.  The  sources  to  be  con- 
sidered will  be  azimuthally  independent.  Let  q and  w 
be  the  displacements  in  the  r and  z directions  and 
define  dilatational  and  rotational  potentials  $ and  y 
such  that 


^r’z)  = H + SrlUr  > 

w(r,z)  = ^ + ^-7  + k>  • 


(c.n 


Let  R2  = z2  r2.  The  P wave  displacement  along  R 


is  then 


Up  = sin0  + wp  cose  = ^ . (C.2) 

We  shall  subsequently  be  concerned  with  the  reduced 
displacement  potential  V,  defined  by 


Up  = 


(5  ■■“•*)  • 


(C.3) 


Therefore , 


-ikaR 


(C.4) 


Point  Compressional  Wave  Source  in  a Half  Space 


Consider  a point  source  of  compressional  waves  in  an 
unlimited  solid.  The  potentials  may  then  be  written 


■ikaR 


0 = A tt 

cw  o R 


(C . 5) 
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Let  such  a source  be  embedded  at  depth  h in  a half 
space.  me  dilatational  potential  is  then 


$ 


cw 


-ikr  (z-h)  'ikr  (z+h) 
e a - P e a 


J (hr) 

o 


dk  , 


a 


(C.6) 


where 


CC.7) 


sin0  = a/c  , 


F 

i 


Jnd  c is  horizontal  phase  velocity  while  6 may  be  identi- 
fied as  the  takeoff  angle  for  waves  emitted  by  the  source. 

The  integral  (C.6)  may  be  evaluated  by  contour  integra- 
tion. The  poles  of  the  integral  give  tie  surface  waves  while 
body  waves  are  obtained  by  evaluating  the  branch  line  contri- 
bution. Restricting  attention  to  the  far  field  (kr  >>  1), 
the  method  of  steepest  descent  may  be  applied.  In  this  case 
the  result  is 


$ = A cos0 

cw  o 


| z-h 


“lkrah  ) -ik(r  z + r) 

p Vttt  j e 


(C.8) 


Consider  the  case  z >>  h which  is  more  or  less  implied 
by  our  far  field  assumption.  Then 


<5 

cw 


= A 


o 


ikr0h 


P e 


-ikr  h 
a 


ik0R 


(C.9) 


I 


i 
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As  h 0 the  coefficient  in  parentheses  approaches 


1 - P 


8 Ve 


i) 


+ 4 Ve 


(C.10) 


A Concentrated  Force  Acting  on  the  Free  Surface  of  an  Elastic 
Half  Space 

Let  us  impose  a vertical  concentrated  force  on  the 
free  surface  in  the  following  way.  Let  the  stresses  at  the 
free  surface  be  given  by 


lTzzl2=o  ' ' If  / J,(kr>  kdk  • 


[T  1 
1 zr J 


z=0 


0 . 


CC.  ii) 


With  this  boundary  condition  the  dilatational  potential 
can  be  shown  to  be  ' 


PF 


oc 


(2k2  - kj)  -ikr  z 

— e J (kr)  kdk 

o 


tut 


(C. 12) 


where  the  Rayleigh  function 


F(k)  = k“[(r2  - i)*  + 4 ror  ] , 


(C . 13) 


and  p is  the  shear  modulus. 

Once  again,  contour  integration  may  be  applied  to  ob- 
tain the  dilatational  potential  for  far  field  P waves: 


PF 


iL 

2TTpk 


r Q 


e 


-ikaR 


CC. 14) 
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where 


Q - 


l)2  + 4 


rare 


(C . 15) 


Relation  Between  a Compress ional  Wave  Source  and  a Point  Force 

It  is  immediately  apparent  that  a compressional  wave 
source  at  vanishingly  small  depth  is  a distinctly  different 
source  type  than  a point  force  acting  downward  on  the  free 
surface.  In  the  first  case  a stress  free  boundary  is  assumed, 
while  in  the  second  a stress  is  imposed  on  the  free  surface. 
For  the  case  of  a whole  space  the  differential  relation  be- 
tween the  two  is  given  below. 

In  an  unlimited  solid  the  dilatational  potential  for 
a vertical  point  force  may  be  written 


4> 


PF 


L 8_ 
4 TT  U) 2 p 92 


-ik  R 
a 

“Tr- 


ee.16) 


where  p is  density. 

Comparing  to  (C.5),  this  implies  that 
i 9$ 


PF 


4irw2pA 


cw 

Sz 


(C. 17) 


That  is,  the  dilatational  potential  for  a concentrated  force 
is  proportional  to  the  space  d .rivative  of  that  for  the  point 
source.  This  result  may  be  generalized  by  taking  the  component 
°f  V4>ci/  in  the  direction  of  the  concentrated  force.  In  fact, 
Eq.  (C.17)  may  be  written 


V2$ 


PF 


4iTpa2A 


n 


V$ 


cw 


(C.'8) 


where  n is  a unit  vector  in  the  concentrated  force  direction. 
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Application  to  an  Explosion  Given  by  a Reduced  Displacement 
Potential  (RDPj  i 


Comparing  Eqs.  (C,4)  and  (C.5),  the  compress ional  wave 

source  amplitude  A is  related  to  the  RDP  by 

0 


A 

o 


= -Y 


CC. 19) 


Therefore,  from  (C.2)  and  (C,9),  the  far  field  displacement 
for  a shallow  source  may  be  written 


U 


cw 


ik  ¥ 
a 


ikr  h 
a 


P 

e 


- ikr  h 
a 


-ik  R 
e a 
R 


(C.20) 


Now  consider  a point  force  at  the  surface.  In  order  to 
write  the  point  force  amplitude  L in  terms  of  a contained 
explosion  quantity  like  the  RDP,  one  must  make  assumptions 
about  the  coupling  of  the  explosion  into  the  ground.  That 
is,  one  must  specify  the  free  surface  boundary  condition  (C.ll) 
in  terms  of  explosion  quantities. 


It  is  reasonable  for  our  purposes  to  assume  that  the 
stresses  on  the  free  surface  are 


(C. 21) 


where  n is  a constant  (<  1.)  and  P is  the  pressure  generated 
by  a contained  explosion  with  RDP  Y,  That  is, 

-ik  R 

M,0  ' " ^ ^ R “ > <C'22) 

where  K is  the  bulk  modulus. 
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From  (C,ll), 


Then 


oo 

M,0  = ’ 


J (kr)  k dk, 
o 


L 9_  e 
3TtF  9z  R 


ik  R 
a 


z = 0 


(C. 23) 


i US  e 
2tt  a R 


-ik  R 
a 


z = 0 


2tt  iv 


cose 


^ nKY 


(C . 24) 


Using  (C.24)  in  (C.14)  and  applying  (C.2),  the  far  field  dis- 
placement is  then 


ik 


U 


a 


PF 


psin  6 


hKQT  | 


-ik  R 
a 


(C . 25) 


Comparison  of  Displacement  for  the  Two  Source  Types 


Now  compare  the  displacement  radiation  patterns  for  the 
point  force  and  the  compressional  wave  source  buried  at  vanish- 
ingly small  depth.  From  (C.20)  the  far  field  displace- 
ment for  a compressional  wave  source  with  h 0 may  be 
written 


U 


cw 


A ¥ 


cw 


-ik  R 
e a 
R 


(C . 26) 


167 


where  Acw  = 1 -P  and  is  given  by  (C.10).  The  corresponding 
displacement  for  the  displacement  due  to  a point  force  at  the 
free  surface  is  given  by  (£.25),  This  may  be  written 


-ik  R 

lka  App  y I 


where 


nKQ 

ysin20 


(C. 27) 


(C. 28) 


For  the  case  of  a2  * 3B2,  (28)  reduces  to 


5 nQ 

^ sin20 


(C.29) 


Comparing  (26)  and  (27) , we  see  that  the  displacements  take 
the  same  form  for  both  sources,  differing  only  by  the  take- 
off angle  dependent  factors  n Anc  and  A . These  factors 

r r cw 

are  tabulated  below  and  are  plotted  in  Fig.  C.l, 


_0 

A 

cw 

App/ n 

0 

0 

0.556 

10 

0.046 

0.554 

20 
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0.549 

30 

0.374 
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40 
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0.536 

50 

0.825 

0.536 

60 

1.000 

0.556 

70 

1.072 

0.627 

80 

0.921 

0.848 

90 

0 

1.667 
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Figure  C.l.  The  takeoff  angle  (6)  dependent  factors  in  the 
far  field  displacement  due  to  a compressional 
wave  source  at  vanishingly  small  burial  depth 
CAcw)  and  a concentrated  force  applied  at  the 
free  surface  (App/p).  The  material  properties 
are  such  that  a2  = 382. 
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Summary 

As  is  seen  in  Fig.  C.l,  the  seismic  waves  directed 
downward  (of  interest  for  teleseismic  monitoring)  are 
distinctly  different  for  a compressional  wave  source  at 
vanishingly  small  depth  and  a concentrated  force  applied 
to  the  surface.  These  two  radiation  patterns  are  at  least 
roughly  representative  of  those  to  be  expected  from  fully 
contained  and  cratering  shots.  In  reality  there  is,  of 
course,  a continuous  limit  between  the  two  as  one  decreases 
the  depth  of  burial  of  a fully  contained  explosion.  How- 
ever, elastically  there  is  no  such  continuous  limit  due  to 
the  different  free  surface  boundary  conditions  for  the  two 
cases . 


* 
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APPENDIX  D 

THEORETICAL  AMPLITUDE  SCALING  FOR  RAYLEIGH  WAVES 
GENERATED  BY  UNDERGROUND  NUCLEAR  EXPLOSIONS 

I.  INTRODUCTION 

An  understanding  of  the  gross  scaling  properties  of 
the  vaiious  kinds  of  waves  generated  by  seismic  events  pro- 
vides a useful  frame  of  reference  for  understanding  the 
actual  data.  As  has  been  pointed  out  by  Bache,  et  al.  (1974), 
the  surface  wave  theory  developed  by  Harkrider  (1964)  predicts 
that  surface  wave  magnitude  should  scale  with  the  product  of 
the  source  material  shear  modulus,  y,  and  the  static  value 
of  the  reduced  displacement  potential,  V . That  is 

Ms  * loS  » (1.1) 

where  is  a consistently  measured  Rayleigh  wave  magnitude. 

The  appearance  of  this  scaling  law  in  the  data  is,  of  course, 
complicated  by  travel  path  and  depth  of  burial  effects,  as 
well  as  by  the  difference  between  the  actual  explosion  source 
and  the  pure  center  of  dilatation  model  _ r which  the  scaling 
law  is  derived.  Still,  our  understanding  of  the  spherically 
symmetric  explosion  contribution  to  the  Rayleigh  wave  is  quite 
helpful  in  separating  out  the  contribution  of  these  other 
effects . 

The  scaling  law  (1.1)  is  not  widely  known  among  those 
interested  in  nuclear  explosion  seismology  and,  indeed,  there 
seems  to  be  a certain  amount  of  confusion  on  this  point.  For 
this  reason  a detailed  derivation  of  the  law  is  given  below. 

We  begin  with  a discussion  of  the  elastic  half  space  solutions 
^ited  by  others  (e.g.,  Marshall  (1970)),  to  provide  a scaling 
law  for  surface  waves.  After  noting  the  logical  inconsistency 
of  this  kind  of  analysis,  we  proceed  to  give  a modified  half 
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space  solution  which  provides  an  elementary  verification  of 
the  scaling  law  (1.1).  Finally,  the  appearance  of  the  same 
scaling  law  for  an  explosion  in  a layered  half  space  is 
briefly  discussed. 

Before  proceeding  with  the  detailed  analysis,  a brief 
discussion  of  the  physical  basis  underlying  this  scaling  law 
is  in  order.  In  terms  of  the  reduced  displacement  potential, 
the  radial  stress  may  be  written 


T 


RR 


!_  + in.  i_  + 

R3  a R2 


(1.2) 


where  a is  the  P wave  velocity,  p is  density,  R is 
the  distance  to  the  point  of  observation  and  the  dot  indi- 
cates differentiation  with  respect  to  the  reduced  time, 
t - R/a.  Then  at  late  time, 


x 


RR  5 


(1.3) 


In  the  same  fashion  it  can  be  shown  that  the  long  time  limiting 
value  of  the  hoop  stress  is 


Tee  15  — ^ 

00  R3  00 


(1.4) 


Therefore,  v/e  see  that  the  stress  in  the  cavity  vicinity 
scales  the  same  way  as  the  teleseismic  Rayleigh  wave.  It 
can  be  said  that  it  is  this  quantity  that  couples  directly 
into  the  Rayleigh  wave. 


172 


* 


R- 2646 


( 


II.  iiLASTIC  HALF  SPACE  SOLUTION 


Consider  a spherically  symmetric  explosion  in  a homo- 
geneous, isotropic  medium.  The  radial  displacement,  U (R,t), 

may  be  written  in  terms  of  the  reduced  displacement  potential 
^(t-R/a)  as  follows: 


UR(R,t) 


a 

' aE 


Ut-R/cQ 

IT 


(2.1) 


where  R is  the  radial  distance  from  the  source,  t is  the 
time  and  a is  compressional  wave  velocity.  Equation  (2.1) 
is  equivalent  to  the  assumption  that  a radial  stress 

tRR  a "P(t)  is  aPPlied  to  the  surface  of  a spherical  cavity 
of  some  radius  a. 


Taking  Fourier  transforms,  (2.1)  may  be  written 


(2.2) 


where  ka  = w/a,  the  compressional  wave  number. 

For  a source  of  this  form  at  a depth  h in  an  elastic 
half  space,  the  surface  displacements  are  given  by  Ewing, 

Jardetskyjmd  Press  (1957),  Eq.  (2.86).  The  vertical  com- 
ponent, cj  is 


wQ(r,w) 


24'  (cj) 


(2k2-k*)  -V  h 
— F(k)  e Jo(kr)  dk  , 


(2.3) 


where  the  Rayleigh  function 


F(k)  = C2k2-k|) 2 - 4k!VjJ  , 
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with 

\>y  = k2  - k*  , y = a or  6 , 

and  kg  * w/3,  the  shear  wave  number. 

The  Rayleigh  wave  contribution  to  the  integral  (2.3) 
is  easily  found  by  evaluating  the  residues  at  the  poles  of 
the  Rayleigh  function.  At  large  epicentral  distances,  the 
Rayleigh  waves  are: 


u 


k?kfi(2kr'kft)  r 

'4i  (2t0  * (<*>)  B exp  -i(v^h+krR-w/4) 


r (kr) 


P'(kT) 


4£(kr) 

31 


(2.4) 


where  kr  = w/cr  is  the  Rayleigh  wave  number  with  cr  the 
velocity  of  Rayleigh  waves.  Also,  v*  = k*  - k*.  There 
appear  to  be  errors  in  the  analogous  expression  of  Marshall 
(1970)  accounting  for  differences  between  the  two.  These 
differences  do  not,  however,  affect  the  scaling  properties. 

The  surface  wave  magnitude  is  measured  from  long 

period  surface  waves.  Therefore,  a scaling  law  is  deduced 
from  the  low  frequency  behavior  of  (2.4).  „ow,  we  note 
that 


lim  iwY(w)  = lim  ¥(t)  = . 

w+0  t-*-°° 

Cherry,  e_t  al , , (1973,  1974)  have  shown  that  iu;?(u))  is 
constant  over  a considerable  frequency  range  and  certainly 
for  those  frequencies  of  interest  for  teleseismic  surface 
waves.  Therefore,  for  low  frequencies 
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0) 


-4  (2tt) 


krk6(2kr-kP 

" Br*5  F'(kr) 


exp 


-i (v'h+krr-n/4) 


(2.6) 


If  the  simplifying  assumption  X = p is  made,  the 
amplitude  of  the  teleseismic  Rayleigh  wave  is 


|w  | s 2.85 
0 


CO  oo 

63/^2 


This  dependence  on  ¥ /B3^2 

oo'  M 

deduced  by  Marshall  (1970). 
law 


(2.7) 


is  essentially  the  same  as  thr c 
It  is  equivalent  to  the  scaling 


M 


s 


st 


log 


¥ 


(2.8) 


which  is  rather  different  from  (1.1). 

The  difficulty  with  the  scaling  law  (2.8)  is  that  the 
analysis  on  which  it  is  based  is  not  relevant  to  the  problem 
at  hand.  The  objective  is  to  determine  the  scaling  of 
Rayleigh  waves  with  the  source  properties,  deleting  travel 
path  effects.  In  the  foregoing,  changes  in  the  source  and 
travel  path  properties  are  linked  and  the  practical  meaning 
of  (2.8)  is  obscure.  In  the  following  section  a modified 

half  space  solution  in  which  this  difficulty  is  surmounted 
will  be  given. 
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We  are  interested  in  the  effect  of  chansin8  source 
region  properties  on  the  Rayleigh  waves  in  a half  space  for 
which  the  material  properties  are  held  constant.  A simple 
moael  providing  interesting  results  is  that  of  a spherical 
inclusion  embedded  in  a whole  space  of  the  half  space  mate- 
rial (Harkrider  and  Bache  (1975)).  The  source  region  geometry 
IS  illustrated  in  Fig.  1.  The  P waves  emanating  from  the 
spherical  inclusion  of  radius  Rs  will  then  be  treated  as  the 
source  for  the  half  space  problem,  with  the  source  region 
assumed  to  be  transparent. 

Referring  to  Fig.  1,  we  begin  by  computing  the  P 
waves  in  material  2 due  to  the  imposition  of  a pressure  pulse 
on  the  spherical  cavity  of  radius  R For  the  spherically 
symmetric  source  in  a whole  space  the  displacement  (UD)  and 
radial  stress  (trr)  written  in  potential  form  are: 

3$. 

\ = air  * 


TRRi  = (Xi+2Mi) 


3 2*i  2\t  34>i 

ap*  R 3R  * 


(3.1) 


where  the  subscript  i - 1,  2,  indicates  the  material  and 

X‘  “ are  the  Lan6  constants.  The  potentials  satisfy  wave 
equations  and  have  the  general  solution  form 


e 'V  eik«  R 

\ R-  + IT  1 - 


* a A n- 

2 2 R 


e^V 

e 2 


(3.2) 
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The  appropriate  boundary  conditions  are  continuity 
of  stress  and  displacement  at  the  material  discontinuity 
and  the  imposition  of  a spectral  pressure  on  the  wall  of 
the  cavity.  That  is, 

tRR[  (RP  ■ • 

Ur  (Rs)  - UR  (Rs)  , (3.3) 

1 2 

trr  C*s>  = trr  Crs>  * 

1 2 

Using  (3.1)  and  (3.2),  the  boundary  conditions  give 
three  equations  in  tha  three  unknown  coefficients  A , B , 
a2-  Since  we  are  interested  in  the  P waves  emanating 1 from 
the  source  region  in  material  2,  A is  of  primary  interest. 
Solving  (3.1)  - (3.3)  this  is 


A 

2 

where 


»i!R)  1 (>  * ‘J.'f  . 
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y e (R) 
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(3.4) 
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with 


y E (R  ) ( n -if7,  (R  ) iF  fR  ) ) 

A = - rp-r1  /b  v ?(u)  *-  e . "He  1 ' 3 / xfF ,(R  (R  )] 
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i 


Using  the  source  specified  by  (3.7)  as  the  transparent 

so jrce  in  a ha.f  space,  the  vertical  surface  displacement  is 
given  by 


<*»  (r,w)  = -2A 

0 2 


D 

/ 


k /2k2 -k2  \ 

J V ~va  h 

— e 2 J (kr)  dk  . (3.8) 


kB  ~ F (k) 
0 2 21- 


Where  r is  epicentral  distance  in  a cylindrical  coordinate 
system. 

Comparing  (3.8)  to  the  analogous  expression  (2.3)  for 
a homogeneous  half  space,  we  see  that  (3.8)  may  be  viewed  as 
the  vertical  surface  displacement  due  to  a source  specified 
y a mod 1 1 led^ reduced  displacement  potential,  y M written 

" I""  °f  ",  and  * «««  transmission  coefficient, 
T(w).  That  is,  set 


with 


Az  (w)  = -y  (w) 


(«)  = T(w)  y (W) 
2 i 


(3.9) 


and 


T(«) 


y E (R  ) (D  "iF.(Rs) 
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The  Rayleigh  wave  contribution  to  the  far  field  vertical 
surface  displacement  is  now  given  by  (2.4)  with  4*  re- 
placing ¥.  2 

The  low  frequency  behavior  of  the  Rayleigh  wave  verti- 
cal surface  displacement  is  needed  to  determine  the  dependence 
of  Mg  on  the  source  parameters.  It  can  be  shown  that 


lim  | T (w 
oj+0 


h*  ■ ycy  - (■  • c y 


(3.10) 


where  the  notation 
3a2 

Q = —L  - 1 
4(32 


has  been  introduced.  If  the  simplifying  assumption  a2  = 3$2 
is  made,  (3.10)  reduces  tu  1 1 


lim  | T (oj)  j 
w-*-0 


s1* 

2 


/ R 


\ 5y 

)+  4 + 


(3.11) 


Then,  combining  (3.11)  with  (3.9)  and  taking  the  ap- 
propriate long  period  limits  (see  (2.5)),  the  amplitude  of 
the  teleseismic  Rayleigh  wave  (analogous  to  (2.7)  for  the 
homogeneous  half  space)  takes  the  form: 
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(3.12) 


where  Tw  is  taken  to  be  the  static  value  of  the  explosion 
reduced  displacement  potential,  ¥ (w) . 

Therefore,  the  dependence  of  on  the  source  para- 

meters for  this  simplified  model  is 

Ms  * log  T YJ  . (3.13) 

This  is  identical  to  the  sc  .ling  law  (1.1),  which  was  our 
objective,  with  the  exception  of  the  factor  T. 

Let  us  now  take  a closer  look  at  what  will  be  termed 
the  ^ transfer  function,  T.  Consider  the  static  problem 
of  applying  a constant  pressure  in  the  cavity  of  radius  R 
in  Fig.  1.  This  problem  can  be  solved  by  taking  the  strain 
potential  of  (3.1)  and  assuming  solutions  in  the  two  mater- 
ials of  the  form: 


IT  * B,  R2 


(3.14) 


Applying  a pressure  of  "P  at  R = R and  requiring  continuity 
of  radial  stress  and  displacement  at  R = Rg , one  can  solve  for 
-he  constant  coefficients.  It  is  easily  shown  that 


A = 
2 


Q l - 


i)  • ( 


1 «■  Q 


( 7 • 1 5) 
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where 


P R3 

o 

4y 


This  result  is  identical  to  the  low  frequency  limiting  case 
for  the  dynamic  problem  given  in  (3.9)  and  (3.10). 

Note  that  in  material  1 the  radial  stress  is  given  by 


trr  (R)  = 4y 

i 


— L + 2QB 
R3  i 


(3.16) 


with  the  uniform  stress  term  being  required  to  maintain 
welded  contact  at  the  material  interface.  It  is  this  term 
which  is  responsible  for  the  presence  of  the  non-unity  factor 
T in  (3.12)  or  (3.13). 

„ From  the  examination  of  the  static  solufon  it  appears 
that  T t 1 is,  at  least  in  part,  a result  of  the  idealized 
geometry  that  has  been  chosen.  The  meaning  of  this  term  for 
the  physical  problem  of  interest  can  be  better  understood  by 
examining  its  value  for  typical  problem  parameters. 

Consider  first  the  case  R#/Rs  s 1.  Then  T s 1 and 

(3.13)  is  identical  to  (1.1).  On  the  other  hand,  when 

Ro/Rs  * 0, 


(3.17) 


Values  of  T for  these  two  limiting  cases  and  two  inter- 
mediate examples  are  plotted  versus  y /y  in  Fig.  2. 
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In  relating  the  simplified  model  studied  here  to  the 
more  general  problem  of  interest,  we  first  note  that  p 
represents  the  elastic  properties  of  the  explosion  source 
region  while  represents  the  "average"  elastic  properties 

of  the  source-receiver  travel  path  for  long  period  Rayleigh 
waves.  The  latter  quantity  is  relatively  constant  over  the 
earth's  surface  and  should  be  approximately  that  for  crustal 
granite,  i.e.,  p^  z 40C-500  kbar.  Nuclear  explosions  have 
been  detonated  in  a variety  of  materials  ranging  from  soft 
alluviums  to  hard  granites.  The  range  of  p /p  values  of 
practical  interest  are  therefore  from  something  less  than 
0.1  to  nearly  1.0.  It  should  also  be  noted  that  a Poisson's 
ratio  of  0.2S  has  been  assumed  for  the  ?;ource  material  for 
Fig.  2.  For  many  earth  materials  Poisson's  ratio  is  some- 
what lower;  for  example,  a value  of  0.20  is  more  appropriate 
for  NTS  tuff.  Assuming  this  value,  the  curves  of  Fig.  2 are 
closer  to  unity.  For  example,  T = 2.0  for  R /R  = p /p  = o 

0 S 1 2 ^ 

and  the  other  values  scale  down  proportionally. 

In  any  case,  for  realistic  explosion  events  it  is 
ve~y  difficult  to  justify  choices  of  the  material  and  geo- 
metric constants  for  our  simplified  model  which  give  values 
of  T much  different  than  unity.  With  this  heuristic  argu- 
ment we  conclude  that  the  appropriate  Rayleigh  wave  scaling 
is 

Ms  = log  [u^J  . (3.18) 
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IV.  THE  RAYLEIGH  WAVES  FROM  AN  EXPLOSION 
IN  A LAYERED  HALF  SPACE 

Harkrider  (1964)  has  given  an  expression  for  the 
Rayleigh  waves  from  an  explosion  at  depth  in  a layered  half 
space.  The  explosion  source  was  modeled  as  a pressure  ap- 
plied to  the  walls  of  a spherical  cavity  of  radius  R just 
as  in  Section  III  above.  That  is,  the  source  potential  is 
given  by  (3.6). 

The  appropriate  expression  for  the  Fourier  time 
transformed  Rayleigh  wave  surface  displacements  due  to  the 
presence  of  this  source  in  a layer  denoted  with  subscript  s 

is  given  by  Harkrider  (1964),  Eq.  (85).  In  our  notation  this 
is 


w.  “ kr  Ks  ar  H<2)(krr)  (4.1) 

where  AR  is  the  amplitude  response  of  the  layered  medium 
and  depends  only  on  the  layered  earth  model  for  the  propaga- 
tion path.  The  remaining  undefined  quantity,  Kg,  is  the 
Rayleigh  wave  excitation  factor  which  depends  weakly  on  the 
source  properties  in  the  following  sense.  The  explosion 
may  be  viewed  as  a discontinuity  in  displacements  and  stresses 
at  the  source  depth.  The  factor  Ks  is  related  to  this  dis- 
continuity and  is  therefore  dependent  on  the  properties  of 
the  source  layer. 

The  Eq.  (4.1)  is  only  rigorously  valid  when  the  mate- 
rial in  the  layered  earth  model  at  the  source  depth  and  the 
material  immediately  around  the  explosion  source  are  the 
same.  This  is  a consequence  of  the  form  of  the  term  K^. 
Therefore,  a direct  derivation  of  the  dependence  of  (4.1)  on 
the  source  parameters  encounters  the  same  difficulty  as  the 
half  space  analysis  of  Section  II;  that  is,  one  is  forced  to 
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change  the  properties  of  the  propagation  path  a long  with  the 
source  material.  There  are  several  way:,  to  circumvent  this 
difficulty:  (1)  Replace  the  material  at  the  source  depth  in 

the  layered  earth  model  with  a thin  layer  of  the  local  source 
material;  (2)  Use  the  simplified  source  region  modeling  of 
Section  III  to  carry  the  explosion  generated  wave  into  the 
average  crustal  model  for  the  travel  path;  (3)  Ignore  the 
inconsistency  and  change  only  V (o>)  in  (4.1)  from  event  to 
event.  All  three  of  these  methods  lead  to  the  same  con- 
clusion; the  scaling  of  Rayleigh  wave  amplitude  with  p ¥ 

A brief  discussion  of  this  result  is  given  below. 


For  the  long  periods  of  interest  for  Mg,  the  ampli- 
tude response,  A^,  and  the  Rayleigh  wave  excitation  factor, 
Kg,  are  essentially  unaffected  by  the  presence  of  a thin 
layer  of  contrasting  material  at  shallow  depths.  This  is 
physically  reasonable  and  has  been  verified  computationally. 
Then  (4.1)  is  equivalent  to  the  scaling  law  (1.1). 

In  Section  III  an  equivalent  explosion  source  was  de- 
veloped for  the  case  of  an  explosion  in  a spherical  region 
of  one  material  embedded  in  a second  material.  This  equiva- 
lent source  can  oe  used  in  the  half  space  solution.  The  re- 
sulting Rayleigh  waves  are  given  by  (4.1)  multiplied  by  the 
source  transfer  function  T(u)  of  (3.9).  The  scaling  law 
is  then  given  by  (3.13).  Just  as  in  Section  III  it  can  be 
argued  tnat  T = 1 for  physically  realistic  problems  and,  once 
again,  the  scaling  law  (4.1)  results. 


The  simplest  procedure  when  dealing  with  the  Rayleigh 
waves  from  explosions  in  different  source  materials  but  for 
which  the  average  properties  of  the  travel  path  remain  con- 
stant is  to  ignore  the  inconsistency  in  the  theory.  The  re- 
sult is  essentially  the  same  as  for  the  more  elaborate  pro- 
cedures outlined  above;  scaling  of  M with  loc  fu  ¥ 1 

S “ 1 
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APPENDIX  E 

NARROW  BAND  FREQUENCY  DISCRIMINANT  CODES 


MAGIC 


MAGIC  (MAGnitude  Indication  Code)  is  flowcharted  in 
Fig.  E.l. 

MAGIC  obtains  its  instructions  from  cards,  and  then 
prints  out  information  to  identify  the  particular  run  being 
performed.  The  code  next  proceeds  to  tape  or  drum,  in 
order  to  input  the  event  time  series  upon  which  the  analysi 
is  to  be  carried  out.  This  time  series  is  usually  a single 
event;  the  code  does,  however,  have  the  capability  of  sum- 
ming several  time  series  in  order  to  artificially  produce  a 
multiple  event  signal. 

If  more  than  one  event  has  been  specified,  MAGIC 
will  add  the  individual  time  series,  each  with  a different 
time  lag  and  scale  factor, 


X(t) 


Z w*  * v . 


in  order  to  produce  a multiple  event  record,  which  is  used 
in  the  rest  of  the  program.  This  system  allows  considerable 
flexibility  in  the  construction  of  synthetic  multiple  event' 
If,  on  the  other  hand,  only  one  individual  even  is  speci- 
fied, then  MAGIC  will  operate  on  this  time  series  alone. 
Where  a "multiple  event"  is  specified  in  the  flow  chart  it 
should  be  kept  in  mind  that  this  could  also  be  one  mdividua 
ivent,  in  which  case  the  first  loop  shown  would  be  inactive. 

The  single  or  multiple-event  time  series  is  then 

narrow-band  filtered.  In  MAGIC,  this  an  be  accomplished 
in  one  of  two  ways. 
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Figure 


Loop  Pack , 
in  order  to 
Follow  this 
Procedure  for 
Several  Pif- 
ferent  Fre- 
quency Filters 


Co  Back,  so 
as  to  do  this 
for  each 
Filter 


t.l.  MAGIC86  Flow;hart  (continued). 
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The  first  filter  method  is  a recursive  method,  imple- 
mented in  subroutine  RFCFIL.  In  this  procedure,  the  time 
series  is  filtered  by  a four-term  recursion  relation  which 
approximates  the  filter 

F(u>)  * jl  ♦ 3-fr2  “jO 

( icju 

in  both  the  forward  and  backward  time  sense,  and  the  two 
results  are  summed,  in  order  to  produce  the  effect  of  a 
phaseless  filter.  This  filter  is  centered  at  frequency  u, 
and  has  a quality  factor  of  Q. 

Another  filter  method  consists  of  performing  a Fast 
Fourier  Transform  (FFT) , multiplying  the  event  in  the  fre- 
quency domain  by  a transfer  function  representing  the  same 
filter,  and  then  inverse  transforming  to  obtain  the  filtered 
signal. 

Regardless  of  which  filtering  method  is  used,  the 
result  will  appear  as  a sinusoidal  carrier  wave  centered  at 
w,  with  an  envelope  containing  information  about  the  event 
time  series. 

MAGIC  then  calls  a subroutine,  PARAMA  which,  following 
a very  complicated  procedure,  locates  the  peaks  of  the 
envelope.  The  peak  of  peaks  at  the  envelope  corresponding 
to  this  particular  w is  also  located. 

The  program  performs  this  filtering  and  maxima- 
selecting  procedure  for  each  of  up  to  20  filters,  and  scores 
the  results. 

MAGIC  next  cycles  through  the  filter  frequencies  m 
again,  and  computes  a spectral  magnitude  for  each.  The 
spectral  magnitude  is  calculated  as 
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M = logio(Ag)  ♦ B (A) 
for  body  waves , and 

M = log  (A  ) + 1.656  log  (A)  + 1.689 
»og  10 

for  surface  waves  , where 

A = event  to  detector  distances  in  degrees, 

B (A)  - distance  correction  function  from  a table, 

A - Ymay  uj/2ttA.  . , 
g MAX  inst’ 

YMAX  ~ maximum  envelope  for  filter  frequency  w, 

Afnst  = amPl^tude  response  of  short-period  LRSM 
seismometer  at  frequency  w . 


MARS 

As  a result  of  experience  with  MAGIC,  a more  rapid  and 
sophisticated  method  of  processing  data  was  devised,  and  this 
formed  the  theoretical  basis  for  the  program  MARS  (Multiple 
Arrival  Recognition  System).  This  code  is  in  the  process  of 
development,  and  those  features  which  are  currently  operational 
are  flowcharted  in  Fig.  E.2. 

MARS  begins  by  reading  instructions  from  cards,  and 
then  reading  a single  time  series  from  drum  or  tape  file. 

No  multiple  event  summation  is  currently  implemented  in  this 
program. 

A FFT  is  performed  on  the  time  series,  and  the 
spectrum  is  displayed.  The  transform  is  preserved  through- 
out the  program,  so  that  only  one  forward  transfsy^  need  be 
done  during  the  run. 
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Figure  E,2,  MARS19  Flowchart 
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MAKS  constructs  the  envelope  oi'  a narrow-hand 
filtered  signal  by  means  of  the  Hilbert  transform  method 
(Bracewell,  1965).  The  Fourier  transform  of  the  input 
signal  is  multiplied  by  the  transfer  function  of  a filter, 
to  give  the  Fourier  transform  of  the  filtered  signal.  This 
in  turn,  is  multiplied  by  -i  sgn(u)  to  produce  the  quadra- 
ture transform.  Both  the  f iltered-signal  transform  and 
the  quadrature  transform  are  then  inverse  FFTed . 

The  filtered  time  series  and  the  quadrature  signal, 
now  both  in  the  time  domain,  are  then  regarded  as  the  real 
and  imaginary  parts  at  a complex  analytical  signal;  the 
modulus  of  the  analytical  signal  is  the  envelope  of  the 
filtered  signal.  An  instantaneous  phase  can  also  be  cal- 
culated for  the  analytical  signal,  and  differentiation  yields 
an  instantaneous  frequency. 

After  these  calculations , the  envelope  peaks  can  be 
located  with  relative  ease.  Each  of  these  peaces  represents 
the  arrival  of  a wave  group  having  the  same  center  fre- 
quency as  the  filter.  All  of  this  can  then  be  repeated  for 
filters  cf  different  center  frequencies.  After  doing  this 
for  a number  of  frequencies,  the  group  arrival  time  versus 
fi  .ter  frequency  can  be  plotted,  in  order  to  obtain  a display 
of  the  arrival  at  dispersed  and  undispersed  wave  modes. 

The  envelopes  produced  by  the  different  filters  can 
be  added  together  in  order  to  observe  the  arrival  of  various 
wave  groups.  If  the  envelopes  are  simply  stacked,  this  will 
tend  to  emphasize  the  arrival  of  undispersed  groups.  On  the 
4 other  hand,  ii  the  envelopes  being  added  are  lagged  relative 

to  one  another  in  the  proper  manner,  then  dispersed  wave 
forms  will  manifest  themselves. 

Since  MARS  does  all  filtering  in  the  frequency  domain, 
♦ this  method  grants  much  greater  flexibility  in  the  choice  of 

filter  shapes  than  was  the  case  in  MAGIC.  Experimentation 
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with  MARS  has  shown  that  the  filter  shape  is  very  important; 
this  point  can  best  be  illustrated  with  pictures. 

Figure  E.3  shows  an  input  time  series,  I.ASA  Bvent  1 507  , 
which  was  obtained  from  ACDA.  This  signal  was  passed  through 
a square  filter,  centered  at  f = 0.5  Hz,  with  quality  factor 
Q ® 10.  (All  of  the  filters  to  be  described  are  shown  in 
Fig.  E.4.J  The  sharp  sides  of  the  filter  are  reflected  in 
the  result.  Fig.  E.5,  where  a pervasive  0.05  Hz  ringing  is 
evident,  which  makes  it  difficult  to  locate  the  arrival  with 
any  certainty. 

When  the  filter  employed  by  MAGIC  is  used,  the  result 
is  as  shown  in  Fig.  E.6.  This  is  a much  more  reasonable 
result,  the  envelope  appearing  only  when  the  signal  does, 
and  with  roughly  the  same  shape.  This  envelope  has  a problem, 
however,  one  which  can  barely  be  seen  in  the  plot:  there 

are  actually  74  maxima  in  this  curve.  The  source  of  this 
problem  is  not  hard  to  find;  considered1  in  the  frequency 
domain,  this  filter  never  actually  goes  to  zero,  no  matter 
how  far  away  from  the  center  frequency  one  looks  at  it.  The 
Event  1507  spectrum  peaks  at  1 Hz,  where  the  amplitude  is 
approximately  6 times  greater  than  it  is  at  the  filter  fre- 
quency of  0.5  Hz.  The  filter  does  not  completely  block  this 
strong  component,  which  is  able  to  leak  through  and  contami- 
nate the  envelope,  with  the  result  described.  The  lesson 
from  this  is  that  the  filter  must  have  a complete  cutoff  at 
some  point  in  frequency. 

When  a semicircul arly-shaped  filter  is  used,  the  re- 
sult is  as  shown  in  Fig.  E.7.  This  looks  quite  a bit  like 
the  square  filter,  and  is  unacceptable  for  the  same  reasons. 

A triangular  filter  brings  about  the  result  shown  in 
Fig.  E.8,  a distinct  improvement  over  the  other  envelopes. 

The  best  results,  however,  come  from  a cusp-shaped 
filter,  the  envelope  from  which  is  displayed  in  Fig.  E.9. 
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Figure  E.3.  LASA  Event  1507,  unfiltered.  Amplitude  (arbitrary 
units)  (take  at  y-axis  numbers). 
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ENVELOPE  FOP  FC  - .500,  o - 10.00000 


TItlE  (SEC) 

Figure  E.6.  LASA  Event  1507,  after  passing  through  MAGIC 
filter  of  fc  = 0.5  Hz,  Q = 10. 
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ENVELOPE  FOR  FC  ■ .500,  0 » 10.00000 


TIME  (SEC) 

Figure  F.7.  LASA  Event  1507,  after  passing  through  semi- 
circular filter  of  f = 0.5  Hz,  Q = 10. 
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Figure  E.8.  LASA  Event  1507,  after  passing  through  triangular 


filter  of  f « 0,5  Hz,  Q = 10. 
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Figure  E.9.  LASA  Event  1507  after  passing  through  cusp. 

shaped  filter  of  f£  = 0.5  Hz,  Q = 10.  P 


I h i s case  possesses  the  optimum  features  of  f.i  ) low 
amplitude  in  the  ;ibsi>ncc  of  signal,  and  ( l>  j immunity  to 
noise  outside  the  passband,  leading  to  a smooth  envelope, 
with  a small  number  of  easily  locatable  maxima.  Of  all 
the  filters  explored  to  date  (there  are  a number  which  have 
not  been  shown  here)  , the  cusp  filter  appears  to  be  the 
best . 
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APPEND l X r 
matched  i- 1 1. j 1 king 

The  COLLAPSE  code,  written  by  D.  B.  Rabenstine  and 
J.  W.  Lambert  of  Teledyne  Geotech  SDL,  was  obtained  and  made 
operational  on  the  S3  Univac  1108.  The  program  is  flow- 
charted in  Fig.  F.l. 

In  the  process  of  adapting  this  code  for  the  S3  com- 
puter, major  restructuring  of  the  main  program  was  performed, 
along  with  numerous  minor  adjustments;  these  alteratiors  do 
not  affect  the  logical  system  described  in  Fig.  F.l,  but  do 
considerably  speed  up  the  execution  of  the  program,  shorten 
the  overall  length  of  Fortran  coding,  and  provide  a small 
amount  of  storage  optimization.  Further  improvement  can  be 
expected  to  include  a probable  speed-up  of  the  Fourier  trans- 
form operation,  v^.ch  forms  the  backbone  of  the  program. 

COL30 , the  latest  version  of  COLLAPSE,  is  capable  of 
drawing  plots  both  on  the  line  printer  and  on  S3's  Gould  500C 
plotter. 

Since  the  COLLAPSE  code,  and  the  matched  filtering 
operation  which  it  performs,  have  been  described  in  detail 
elsewhere  (Alexander,  et^  aj^,  1971),  the  remainder  of  this 
discussion  will  be  concerned  only  with  a demonstration  of 
the  program's  ability. 

Figure  F.2a  shows  Event  6545,  a presumed  explosion  in 
eastern  Kazakhstan,  as  recorded  at  the  Oyer  subarray  in 
Norway.  A two  second  time  window,  commencing  with  the  start 
of  the  P-wave  signal  at  approximately  t = 20  seconds,  was 
extracted  and  used  as  a filter.  Figure  F.2b  is  for  Event 
6550,  also  a presumed  Kazakhstan  explosion,  which  is  sub- 
merged in  strong  microseismic  noise;  the  explosion  itself  is 
indicated  by  an  arrow. 
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(a)  TIME  (SEC)  CIO  5 IDJ  TIME  tSEC)  CIO  5 

CONVOLUTION  OF  FILTER  VS4S  WITH  SIGNAL  4SS0 


TINE  CSEC> 


Figure  F.2.  Matched  filter  test  for  two  presumed  Russian  ex- 
plosions. A portion  of  the  event  in  (a)  was  used 
as  a filter  for  the  event  in  (b) , resulting  in  the 
response  shown  in  (c) . 
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When  the  mate  lied  filtering  operation  is  performed,  l lie 
result  (Fig.  F.2c)  shows  a very  clear  and  significant  response 
to  the  hidden  event,  at  approximately  t — 17  seconls. 

Figure  F.3a  is  the  time  series  of  an  Asian  earthquake, 
Event  6224.  When  this  is  used  as  a filter  for  Event  6550 
(Fig.  F.3b),  the  result  (Fig.  F.3c)  indicates  no  significant 
response . 

This  demonstrates  the  ability  of  the  matched  filtering 
operation  to  discriminate  between  different  classes  of  events. 
The  principal  question  which  remains  concerning  matched  filter' 
ing  consists  of  estimating  the  confidence  level  (or  false 
alarm  probability)  associated  with  a given  level  of  response. 
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(a)  time  csec)  cio  ) (b)  time  cseouo  ) 

CONVOLUTION  OF  FILTER  <22<  nITH  SIGNAL  '<550 


(C)  TME  tSECl 


Figure  F.5.  Matched  filter  test  for  an  apparent  Russian  ex- 
plosion (b),  filtered  by  an  Asian  earthquake  (a) 
yielding  the  result  shown  in  (c) . 
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illl  l.(vMI  I VAl.liNT  lil.ASTIC  SOURCE  Dl.TI.KM J NAT  1 ON  I OK 
NONI.INI.AU,  FINITE  DIFFERENCE  EARTHQUAKE 
AND  EXPLOSION  SOURCE  CALCULATIONS 

G.l  INTRODUCTION 

Powerful  nonlinear  finite  difference  computational 
techniques  can  be  applied  to  obtain  a realistic  deterministic 
representation  of  an  earthquake  or  explosion  source  of  seismic 
energy.  In  order  to  link  these  nonlinear  numerical  calcula- 
tions with  efficient  elastic  wave  propagation  techniques,  it 
is  necessary  to  represent  the  souice  in  terms  of  an  equiva- 
lent elastic  source.  Such  a representation  in  terms  of  an 
expansion  in  spherical  harmonics  was  first  suggested  for 
geophysical  applications  by  Archambeau  (1964,  1968).  For 
any  arbitrary  nonlinear  source  calculation,  this  equivalent 
elastic  source  representation  can  be  derived,  if  the  source 
calculation  can  be  carried  into  the  small  displacement, 
elastic  regime. 

The  linkage  between  nonlinear  numerical  explosion 
calculations  and  elastic  wave  propagation  techniques  via  an 
equivalent  elastic  source  has  been  routinely  accomplished  at 
S3  for  spherically  symmetric  explosions  (e.g.,  Cherry,  et  al., 
(1973,  1974b)).  A much  more  complex  source  configuration,  an 
explosion  in  an  axisymmetric  tunnel,  has  also  been  treated 
(Cherry,  e^t  al . , (1974a))  with  the  equivalent  elastic  source 
determined  from  the  output  of  a 2-D  axisymmetric  Lagrangian 
finite  difference  code  which  had  previously  been  linked  to 
an  Eulerian  hydrodynamic  calculation  of  the  explosion. 
Recently,  a 3-D  deterministic  simulation  of  a stick-slip 
earthquake  faulting  has  been  incorporated  into  a Lagrangian 
finite  difference  code  (Cherry  (1974)).  In  order  to  obtain 
an  equivalent  elastic  source  for  this  model,  it  was  necessary 
to  construct  a code  of  complete  generality,  that  is,  with  no 
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of  this  appendix. 

We  have  developed  a series  of  general  purpose  computer 
programs  that  determine  the  equivalent  elastic  source  (EES) 
representation  for  any  arbitrary  nonlinear  numerical  source. 
The  MULTEES  (M'JLTipole  Coefficients  for  an  Equivalent  Elastic 
Source)  scries  of  programs  calculates  the  multipole  coeffi- 
cients for  any  general  numerical  (earthquake  or  explosion) 
source  for  which  the  calculation  has  been  carried  into  the 
small  displacement  elastic  region.  MULTEES  has  been  designed 
to  process  the  output  generated  by  either  two  or  three 
dimensional  source  codes.  Tims,  the  same  programs  can  be 
applied  to  underground  explosions  and  earthquakes.  Further- 
more, the  numerical  source  codes  linking  to  MULTEES  can  be 
configured  for  any  arbitrary  source.  No  special  symmetries 
are  demanded. 

The  primary  mathematical,  numerical  analysis,  and 
computational  techniques  underlying  MULTEES  will  be  described 
in  the  following  sections  along  with  a narrative  outline  of 
how  MULTEES  operates. 
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1 he  representation  of  a spatially  limited  source  in 
tern’s  of  a multipole  expansion  in  spherical  harmonics  has 
been  discussed  in  a number  of  previous  S3  reports  (e.g., 
Cherry,  et^  al . (1974a),  Bache , et_  al . (1974)).  However,  in 
order  to  introduce  the  definitions  and  notation  required  for 
describing  MULTEES,  the  derivation  of  the  expansion  (multi- 
pole) coefficients  in  terms  of  potentials  will  be  reviewed 
here . 

The  equation  of  motion  in  a homogeneous,  isotropic, 
linearly  elastic  medium  may  bo  expressed  as 


— - = v*  7 (7*u)  - v27  X (7  X u)  , 
3t2  S 


(G.l) 


where  u(r,t)  is  the  particle  displacement,  while  vp  and 
vs  are  the  compressional  and  shear  wave  velocities.  Define 
the  scalar  and  vector  potentials  by 

X = 7 • u , 

4 ** 


X = j V x u . 


(G.2) 


By  introducing  these  potentials,  the  equation  of  motion  re- 
duces to 


7 " vp  7X||  - 2v2  7 x x • 


(G.5) 


Taking  first  the  divergence,  then  the  curl  of  Eq.  (G.3),  the 
scalar  and  vector  potentials  may  be  shown  to  satisfy  the  wave 
equation  for  an  elastic  medium, 


(v2  - — ~ Xa  (r,t)  = 0 , 
\ v2  3t2/  a ** 


(G.4) 


211 


R- 2646 


where  u - l,2,S,4  w i t li  v v.  Ini  | I , , S .uul  v v . 

.IS'  i,  1 1 

To  derive  solutions  to  Lq.  (li.4),  introduce  the  Fourier 
transform  of  the  potentials  defined  as 

00 

X0(r,w)  = / Xa(£,t)  e"la)t  dt  ; (G.S) 

.00 

conversely,  by  means  of  the  Fourier  integral  theorem,  x can 
be  expressed  in  terms  of  x as 


Xa(r,t) 


2 it  f Xa(r.“) 

- 00 


(G.6) 


For  convenience  herein,  let  F ( f ) represent  the  Fourier 
transform  of  the  function  f;  that  is, 


00 

F [f (t) ] = f(w)  » f f(t)  e'iwt  dt  . 

- 00 


(G.7) 


Thus,  Eq.  (G.4)  reduces  to  the  wave  equation  for  the  Fourier 
transform,  namely, 

(v2  + kj)  x0(£,«)  = 0 , (G.S) 

where  the  wave  number  ka  = w/va«  By  this  means,  the  wave 
equation  for  elastic  body  waves  (G.4)  has  been  reduced  to 
the  standard  Helmholtz  equation  (G. 8),  which  is  separable  in 
spherical  coordinates.  The  ’’outgoing  wave"  solution  for  xa 
can  therefore  be  expanded  in  terms  of  spherical  eigensolutions 
of  Eq . (G . 8)  : 
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Xa(r,6,0,w)  = 


a i 

E El 

£ = 0 m=0  s=0 


hl2)(kar)  *!:!(“)  V 


im s(u)  («.9) 


where  A<“>  are  multipole  coefficients,  h<2>  are  spherical 
Hankel  functions  of  the  second  kind,  and  are  (un- 

normalized)  spherical  harmonics  (Morse  and  Feshbach  (1953)). 
In  terms  of  the  notation  used  in  earlier  S3  reports  e.g. 
Cherry,  et  al.  (1974a),  A<“>=aC“)  and  A$a> 

For  these  discussions  the  following  definitions  will  sufrice: 


,(2) 


(kar)  = J£(kar)  - iyA(kar)  , 


v'G.10) 


,m 


Jims 


(e  ,❖) 


P^cosG)  cosm$,  s = 0 


P™(cos6)  sinm<}>,  s = l 


CG. 11) 


where  j % and  y£  are  spherical  Bessel  functions  of  the 
first  and  second  kind,  respectively,  and  P™(cos0)  are 
associated  Legendre  functions.  The  orthogonality  of  the 
(unnormalized)  spherical  harmonics  is  expressed  by 

2 tt  IT 

/ / YJlms(e,<<))  sin6  de  d<> 

o •'o 

Km  6l'l  6m"m  6s's  * °) 

= Ka  VjI  V0  6s%  6s0  (m  = 0) 

where  the  normalization  factor  N is 

Jim 


CG. 12) 


N 


(47;/em)  (A+m) ! 


4ra  = TH  + l)  (Jl-m) ! 


(G . 13) 


v ith  cq  = l and  cm  = 2,  m ^ 0 
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Xa(r,6, «!>,«)  = 


E E E «{2)Ik,D  *, 


(6 ,4>)  (G.9) 


JL  = 0 m=0  s = 0 


*v  r cx  i r?'! 

where  Akm'  are  multipole  coefficients,  J are  spherical 

Hankel  functions  of  the  second  kind,  and  Y are  (un- 
normalized) spherical  harmonics  (Morse  and  Feshbach  (1953)). 
In  terms  of  the  notation  used  in  earlier  S3  reports,  e.g., 
Cherry,  et  fil*  (1974a),  * a£>  and  ijjj  »«. 

For  these  discussions  the  following  definitions  will  suffice: 


hl2)CV)  ’ Mkar)  ' • 


(G. 10) 


(6,0) 


P™(cos0)  cosm0,  s = 0 
P™(cos0)  sinm0,  s = 1 


(G. 11) 


where  j ^ and  y ^ are  spherical  Bessel  functions  of  the 
first  and  second  kind,  respectively,  and  P™(cos0)  are 
associated  Legendre  functions.  The  orthogonality  of  the 
(unnormalized)  spherical  harmonics  is  expressed  by 


2tt  tt 


ff 

n •'r\ 


Y£'m's^6  Y)lms  sin0  de  d<(> 


0 -'o 


= N?  60,0  6 , 6 . (m  ? 0) 

JlmJlJimmssv'  1 

33  N*  6 * 6-  6 (m  = 0) 

J6  o mo  So  So  ' 


(G. 12) 


where  the  normalization  factor  N.„  is 

Jim 


(4tt/c  ) (Jl+m) ! 

v2  _ nr*  ' 

Jim  fMI  U-m) ! 


(G . 13) 


with  e = 1 and  e 3 2,  m / 0. 
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Application  of  the  orthogonality  of  the  spherical 
harmonics  to  Eq.  (G.9)  yields  the  following  expression  for 
the  multipole  coefficients: 


c-) 


2 7T  7T 


Ntm  ^2)CkaR) 


f f Xa(R»9,$,w) 

Jc\ 


r0  •'O 


x ^ms^6*^  sin6  de  d(J>  • 

Substituting  Eq.  (G.5)  for  jj  in  the  above 
that 


CG. 14) 


expression  shows 


:[A£2<R- 


t) 


hPV0R) 


00 

/ a£s(r. 


t)  e-lut  dt 


CG. 15) 


where 


2tt  tt 


‘ ZT  f f xa(R,e,*,t) 

Jtm  J0  •/0 


Jims 


(6,<J>) 


x sinQ  d6  d<J>  . 


(G . 16) 


The  quantities  AW(R,t)  defined  b>  Eq.  (G.16)  are  referred 
to  as  the  multipole  coefficients  in  the  time  domain.  An  in- 
terpretation of  can  be  derived  in  the  following 

manner.  Starting  with  the  expression  for  xa  given  by 
Eq.  (G.6),  replace  xa  by  its  multipole  expansion  Eq.  (G.9) 
Then  substitute  (G.15)  for  jW.  Applying  the  Fourier 
integral  theorem  then  yields  in  the  time  domain, 
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(C.I7) 


&=0  m=0  s=0 


which  is  analogous  to  the  muWpole  expansion  of  x (r,w) 
given  in  (G. 9) . 

'•3  THE  numerical  CALCULATION  OF  THE  MULTIPOLE  COEFFICIENTS 
TMULTEESl  " “ — — 


The  numerical  calculation  of  the  mul.-.ipole  coefficients 
for  the  equivalent  elastic  source  centers  on  the  following 
set  of  integrals: 


A O) 

Jims v J 


00  / 2tt  7T 

nTTr  J li^~  I f Xc(R>0>*-t)  sin0  de  d * 

1 a J -°°(^Jlm  0 JQ 


x e dt  . 


CG.18) 


The  form  of  these  integrals  shows  that  the  equivalent  elastic 
source  representation  for  any  arbitrary  nonlinear  earthquake 
or  explosion  source  can  be  derived,  if  the  numerical  calcula- 
tion can  be  carried  to  some  radius  R in  the  elastic  region. 

In  MULTEES,  the  integral  (G . 18)  is  carried  out  in  two 

steps : 

1.  Double  numerical  integration  on  the  surface 
of  a sphere  of  radius  R for  each  time  at 
which  the  potentials  are  sampled.  This  spatial 
or  area  integration  yields  the  multipole  co- 
efficients in  the  time  domarn,  A^fR.t)  as 
given  by  Eq.  (G. 16) . 

2.  Integration  over  time  to  obtain  the  coeffi- 
Cients 
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cumput.i  t iou.i  1 problem  Mini  is  carried  out  using  lust  louiier 
trails  l orm  techniques. 

Computational  step  (1),  the  calculation  of  the  multi- 
pole coefficients  in  the  time  domain  defined  by  Eq.  (G.16), 
requires  the  evaluation  of  a set  of  double  integrals  over 
the  surface  of  a sphere  of  radius  R.  These  double  integrals 
are  of  the  general  form 


A(R) 


/ J 


F(R,0,<}>)  sine  d6  d$  , 


where  p(R,6,<f>)  = 


^t2  ^ Jims  ^ ® ^ ) • While  the 

■elBell  behaved  fnnrti  on  c rn  n 


YJlms^0,<})^NJm  are  ^e11  behaved  functions,  the  xa(R>9,<fO 
are  to  be  obtained  from  discreetly  sampled  values^n  a 
finite  difference  grid.  In  general,  these  values  will  be 
available  at  a limited  number  of  non-uniformly  spaced  points 
in  the  neighborhood  of  the  spherical  surface  of  radius  R. 

The  program  MULTEES  evaluates  the  integral  (G.19)  by  first 
performing  a spatial  interpolation  to  obtain  an  evenly  spaced 
discreet  representation  of  F(R,0,$)  on  the  spherical  sur- 
face, then  performing  a numerical  double  integration.  Before 
describing  the  computational  techniques,  it  is  necessary  co 
discuss  briefly  the  form  of  the  finite  difference  code  out- 
put expected  bv  MULTEES. 
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Cl . 4 MULTI.  i.S  INPUT 

MULTKES  is  designed  to  operate  with  any  numerical 
source  code  employing  a two-  or  three-dimensional  rectangu- 
lar grid  system.  For  compatibility  with  codes  employing 
spherical  or  cylindrical  coordinate  systems,  only  minor 
modifications  are  required.  In  fact,  the  majority  of  the 
special  techniques  formulated  in  developing  MULTEES  were 
required  for  compatibility  ^ith  rectangular  grid  systems. 

Inasmuch  as  the  multipole  coefficients  are  defined 
in  terms  of  potentials,  it  is  necessary  that  the  source 
calculation  generate  either  values  of  xa  directly  or  re- 
lated variables  from  which  xa  can  be  derived.  The  poten- 
tials were  defined  in  terms  of  the  divergence  and  curl  of 
the  displacement  field  by  (G.2).  In  its  current  form, 

MULTEES  accepts  the  divergence  (x  ) and  0.5  times  the  curl 

4 

(x)  as  generated  output  fiom  the  source  calculation  at 

•v 

specific  points  of  interest  in  the  elastic  region. 

The  determination  of  these  selected  points  of  interest 
will  be  described  below,  but  it  is  first  necessary  to  dis- 
cuss another  point  that  should  be  addressed  :n  order  to  en- 
sure the  valid  linking  of  MULTEES  to  a source  calculation. 

The  discussions  below  will  be  carried  out  with  general 
reference  to  thiee  spatial  dimensions,  since  the  primary 
objective  in  developing  MULTEES  was  to  make  available  a 
general  purpose  code  for  computing  multipole  coefficients 
for  any  3-D  numerical  source. 

In  a finite  difference  code,  saved  and  computed  values 
are  computed  either  at  node  points  or  at  zone  centers.  The 
potentials  are  zone-centered  quantities  and  therefore  are 
available  at  the  centers  of  the  rectangular  prisms  or  cells 
formed  by  the  grid  system.  This  fact  is  important,  first 
for  determining  the  points  at  which  the  potentials  are  to 
be  saved  by  the  source  calculation;  and,  secondly,  in 
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communicating  to  the  .source  code  the  idem  i T i cat  ion  of  the 
stations  at  which  these  variables  are  to  be  monitored. 

In  order  to  maintain  the  proper  spatial  definition  of 
source  variables,  such  as  displacements,  velocities  and 
potentials,  the  rectangular  grid  or  "source  coordinate  system" 
on  which  MULTEES  is  based  must  be  identical  with  that  used  by 
t e numerical  source  code.  This  is  accomplished  in  MULTEES 
by  applying  the  same  algorithm  to  generate  the  rectangular 
grid  system  that  was  used  in  the  source  code  iteself.  The 
grid  system,  which  need  not  be  evenly  spaced,  is  defined  by 
three  one-dimensional  real  arrays:  CSX  (NX),  CSY(NY)  and 

CSZ(NZ).  A grid  pci,,*  i5  specified  by  a s’;  of  'J; 

indexes  (I.J.K)  where  the  X.  Y,  and  Z coordinates  of  the  point 
are  defined  by  CSX(I),  CSY(J),  and  CSZCK),  respectively.  I„ 
order  to  apply  the  general  algorithm  referred  to  above  the 

following  five  input  specifications  must  be  given  for  all 
three  coordinates: 

CSX(l)  - location  of  the  first  point  along  the  X-axis; 


DELX 


SCALEX 


NSCALEX 


total  number  of  grid  points  on  the  X-axis; 

the  original  or  constant  interval  spacing 
between  grid  points  on  the  X-axis; 

= a scale  factor  that  is,  multiplied  by  the 
previous  interval  spacing  in  order  to 
introduce  "stretching"  of  successive  grid 
points  on  the  X-axis;  for  example,  if 
SCALEX  = 1.10,  and  DELX  = 1.0,  the  spacing 
between  successive  grid  points  would  be 
1.000,  1.100,  1.210,  1.331,  1.464,  etc.; 

the  grid  point  index  that  identifies  the 
start  of  "stretching"  of  the  grid  points 
on  the  X-axis. 
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in  advance  the  set  of  stations  to  be  monitored  during  the 
numerical  source  calculation.  A ’’monitoring"  or  "save" 
station  is  identified  by  a set  of  three  indexes  (I,J,K)  that 
define  the  point  in  the  rectangular  grid  at  which  the  X-, 

Y- , and  2-components  of  the  displacement  field  and  velocity 
are  to  be  saved.  In  addition  to  the  displacement  and  velocity 
components  at  the  save  station  grid  point,  the  source  calcula- 
tion will  also  sample  the  four  potentials  xa  (a  = 1,2, 3, 4) 
defined  at  the  centers  of  the  eight  cells  surrounding  the 
point  (I,J,K).  Following  the  labeling  convention  that  any 
given  cell  is  identified  by  the  three  indexes  of  the  grid 
point  with  the  three  highest  integer  values,  the  station 
(I » J , K)  will  sample  Xa's  at  the  centers  of  the  following 
eight  cells: 

1.  1 + 1,  J+l,  K+l 

2.  I,  J+l,  K+l 

3.  I,  J,  K+l 

4.  1+1,  J.  K+l 

5.  1+1,  J+l,  K 

6.  I,  J+l,  K 

7.  I,  J,  K 

8.  1+1,  J,  K . 

In  short,  the  key  information  that  MULTEES  provides 
to  the  source  code  is  ? list  of  monitoring  stations.  The 
corresponding  output  that  the  source  code  then  generates 
for  MULTEES  during  each  time  cycle  is  the  set  of  saved 
variables  for  each  monitoring  station.  There  are  a total 
of  38  saved  variables  for  each  station  per  time  cycle 
generated  by  the  source  calculation,  namely, 
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1.  X-,  Y-,  aiul  --components  of  the  d i sp  1 a coition  I 

field  (u,  v,  and  w) ; 

2*  X-,  Y-,  and  Z-components  of  the  velocity 

(u , v,  and  w) ; 

3*  xa  at  the  center  of  the  8 sur- 

rounding cells  in  the  sequence  defined  by  the 
convention  given  above. 

In  addition,  the  cycle  number  and  total  elapsed  time  are 
also  saved  with  the  above  variables  for  each  time  cycle  the 
source  calculation  is  run. 


G.5  DOUBLE  NUMERICAL  INTEGRATION  METHOD 

By  a double  application  of  the  familiar  one-dimensional 
integration  formula,  known  as  Simpson's  rule,  the  integral 


XB  YB 


f f 

aA  yYA 


F(X,Y)  dX  dY 


CG . 20) 


may  be  approximated  by 


where 


(hk/9)  [(F  + F + F + F ) 

13  31  33 

+ 4 (F  + F +F  + F ) + 16  F 1 . 

12  21  23  32J  22j 


(G . 21) 


X(I)  = XA  + (I-l)h 
Y (J)  = YA  + (J-l)k 
Fjj  = F(X(I) , Y(J))  . 

The  numerical  integration  of  the  double  integral  Eq.  (G.20) 
has  thus  been  reduced  to  a simple  two-way  integration  over  a 
rectangle  using  only  two  evenly-spaced  intervals  in  both  X 
and  Y. 
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H>'  apply ii»K  till'  formula  given  in  l.q.  (G.21)  to  sul»- 
rcctangles , and  adding  the  results,  the  two-dimensional 
generalization  of  the  parabolic  rule  can  be  derived  in  the 
form,  as  given  in  Hildebrand  (1956) 


XB  YB 

f f F(X,Y)  dX  dY 


'XA  YA 


2 

(HX  • 

HY/9)  [(F 

+ 4 F +2 

F 

i l 

2 1 

3 1 

+ 

4(F  , 

+ 4 F +2 

F + . . . + 

F 

1 2 

2 2 

3 1 

n2 

+ 

2 (F 

+ 4 F +2 

F + . . . + 

F 

1 3 

2 3 

3 3 

n3 

+ 

...  + 

(F  + 4 F 

+ 2 F + . 

im 

3m 

• • 

F ) 

ni } 


+ F ) ] + E , 
nm'J  » 


(G. 22) 


where 


HX*HY 


n (HX) 


a^cr,^) 

3X4 


m(HY) 


34rcr  ,n  ) 
**  2 2 


3Y* 


HX  « (XB-XA) / (NX-1 ) 

HY  - (YB-YA)/ (NY-1) 

FU  5 F C1 , - F (X  (I)  , Y (J)  ) 

NX  = n - number  of  points  in  the  one-dimensional 
X integration 

NY  = m = number  of  points  in  the  one -dimensional 
Y integration 

XB,  XA  = upper  and  lower  limits  of  integration  for 
the  X variable 

YB,  YA  = upper  and  lower  limits  of  integration  for 
the  Y variable 

lie  in  the  interval  (XA,XB),  and 
r,j  » h2  lie  in  the  interval  (YA, YB) . 
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Both  NX  and  NY  are  odd  integers,  not  necessarily  equal 
corresponding  to  a division  of  the  rectangle  XA-Xll-YB-YA 
into  an  even  number  of  subrectangles.  A major  consequence 
of  using  evenly  spaced  intervals  (HX  and  HY)  in  the  double 
integration  formula  is  the  reduction  of  the  error  term  to 
the  analytic  expression  given  above.  A corresponding  formula 

can  be  derived  for  non-equal  spacing,  but  the  error  will  not 
be  so  readily  known. 

In  order  to  perform  the  numerical  evaluation  of  the 
required  double  integrals  defined  in  Eq . (G.16),  MULTEES 

applies  the  formula  Eq.  (G.22)  to  approximate  these  integrals 
by  multiple  sums  of  the  following  form: 


( 
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F ( I » J) 


F (I » J) 


F (I , J) 


F(I,J)  is  not  an  analytic  function  but  represents  values  of 

the  integrand  at  the  point  (I,J)  that  must  be  available  in 
tabular  form. 

Before  this  double  numerical  integration  formula  can 
be  applied,  the  "integration  points"  (I,J)  on  the  surface  of 
the  sphere  of  radius  R must  first  be  determined.  MULTEES 
offers  the  user  the  option  of  specifying  both  the  integration 
limits  and  the  number  of  integration  points  to  be  used  by 
means  of  the  input  variables: 
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0 >* nd  <p  integrations,  respectively. 


°i’  92’  *x»  ^2  = lower  and  upper  limits  in  the  6 

and  <p  integrations,  respectively. 

MULTHES  calculates  the  interval  spacings  and  the  integration 
points  (in  radians)  by  means  of  the  definitions  given  above 
for  HX,  HY,  X(I),  and  Y(J).  The  integration  points  for  the 

6 and  4>  integrations  are  stored  in  the  one-dimensional 
real  arrays 


G 


r 


♦ ( 


THETA(Nl)  and  PHI(N2), 

so  that  the  integration  point  (I , J)  is  given  by 
THETA(I) , PHI(J). 


The  total  number  of  integration  points  is  the  simple  product 
of  N1  and  N2.  Both  N1  and  N2  must  be  odd  integers,  though 
not  necessarily  equal.  To  assist  the  user  in  selecting 
values  for  N1  and  N2,  MULTEES  also  calculates  the  number  of 
points  'odd  integer)  in  both  the  0 and  4>  integrations 
corresponding  to  an  angular  interval  having  an  arc  length 
approximately  equal  to  the  linear  spacing  between  successive 
grid  points  in  the  rectangular  grid  system.  In  other  words, 
the  largest  odd  integer  is  computed  corresponding  to 

R • (9  - 9 ) /DEL  , 

2 1 


where  R is  the  radius  of  the  sphere,  6 and  6 are  the 
upper  and  lower  limits  of  the  integration  over  (^and  DEL 
is  an  average  linear  spacing  between  grid  points  used  in 
the  source  calculation. 


In  order  to  apply  MULTEES  in  a two-dimensional  case, 
one  of  the  input  variables,  for  example,  N2,  must  be  set 
equal  to  1.  The  effect  of  setting  N2  = 1 is  to  reduce  the 
double  (area)  integration  over  a spherical  surface  to  a one- 
dimensional  or  single  integration  along  an  arc  on  the  sphere. 


I 


..  . .. 
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ilio  double  i ii  t eg  r;i  I ion  routine  in  MiJI.TI  IS  i r.  Jt-‘.  a ^ n«*  J i „ 
handle  automatically  such  an  event  and  to  produce  the  correct 
numerical  value  of  the  corresponding  single  integral.  The 
only  adjustment  that  is  made  by  MULTEES  when  N2  = 1,  is  the 
removal  from  Eq.  (G.23)  of  the  factor  2 H2/3,  where  H2  is 
the  interval  spacing  for  the  4 integration.  In  other  words 
when  N2  = 1 the  double  integration  routine  yields  (2/3)H2 
times  the  result  of  applying  the  one -dimensional  parabolic 
rule  for  the  0 integration. 

MULTEES  is  designed  to  perform  the  required  double 
numerical  integrations  on  the  sphere  between  any  arbitrary 
set  of  angular  limits  - 0^  and  4j  - ^ as  specified  by 
the  user.  The  purpose  for  offering  this  flexibility  is  to 
allow  the  user  to  set  up  the  operation  of  MULTEFS  in  the  most 
efficient  configuration  that  is  suitable  for  the  conditions 
pertinent  to  the  problem  at  hand.  If  the  source  displacement 
field  possesses  certain  symmetry  properties,  one  or  both  of 
the  one-dimensional  angular  integrations  can  he  reduced  to 
a shorter  interval.  In  such  cases,  the  r.ource  calculation 
itself  would  be  restricted  to  a fraction  of  full  space,  per- 
haps only  a hemisphere  or  a single  oct, nt  of  a sphere.  In 
the  latter  case,  the  limits  on  both  0 and  4 might  there- 
fore be  reduced  correspondingly  to  0.0  and  tt/ 2 . While  the 
default  values  in  MULTEES  are  0 to  tt  for  0 and  0 to  2n 
for  4,  the  system  offers  the  user  convenient  options  cor- 
responding to  other  frequently  requested  configurations. 
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In  order  to  derive  values  of  the  integrand  in  Lq.  ((1.18) 
at  the  prescribed  integration  points  on  the  surface  of  the 
sphere,  given  a set  of  values  in  the  neighborhood  of  the 
sphere,  a three-dimensional  ini.erpolation  method  is  applied. 

If  the  values  of  a function  f(Xp,Yp,Zp)  are  available  at  the 
eight  nodes  defining  a rectangular  prism,  then  the  value  of 
f(X,Y,Z)  at  any  spatial  point  within  the  prism  can  be  derived 
according  to  the  following  interpolation  formula  (see 
Zienkiewicz,  Finite  blement  Methods  in  Engineering  Science. 

P.  121), 


8 

- £ f«p,n  ,t  )S(P)  , (G-24) 

P-l 

where  the  shape  functions  for  a "linear"  element  (8  nodes) 
rectangular  prism  are  given  as 


S(p)  = (1/8) (1  + SSp)(l  ♦ nnp)(l  + CCp)  . (G.25) 


The  set  (£»n,C)  specify  a point  in  a normalized  coordinate 
system  for  the  rectangular  prism,  defined  by  the  following 
eight  nodes: 


Node 

1 

2 

3 

4 

5 

6 

7 

8 


1 

-1 

-1 

1 

1 

-1 

-1 

1 


1 

1 

-1 

-1 

1 

1 

-1 

-1 


1 

1 

1 

1 

-1 

-1 

-1 

-1 


CG. 26) 
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^°r  each  integration  point  (X.Y.Z),  a transformation  must 
first  be  made  to  the  normalized  coordinate  system  (C,n,C) 
according  to  the  following  prescription; 


2 (X  - X ) 

F = L_ 

^ X - X 
1 2 

2 (Y  - Y ) 


W 


2(Z  - Z ) 
c - *-*- 


(G. 27) 


where  the  origin  of  the  rectangular  prism  is  given  by 


Xo  “ Xi  ' ' V/2 

Yo  “ \ " - Y„)/2 

Zo  = Z,  - <Z,  - Zs3/2 


(G. 28) 


and  the  subscripts  specify  the  nodes  of  the  prism  (P  = 1,2, 

...  8),  In  other  words,  the  first  node  of  the  prism  is  given 
by  (X^Y^Z^)  in  rectangular  coordinates  and  by  (1,1,1)  in 
the  normalized  coordinate  system  having  an  origin  at 

(Xa’Yn’ZJ. 

0 0 0 

To  illustrate  how  the  interpolation  formula  is  applied 
in  MULTEES,  the  required  computational  operations  performed 
during  pre-processing  will  be  outlined  below.  Since  the 
source  calculation  is  performed  in  a rectnagular  grid,  each 
integration  point  (I , J)  defined  by  THETA (I)  and  PHI  (J)  must 
first  be  transformed  to  rectangular  coordinates  (X,Y,Z).  From 
the  coordinates  of  the  grid  system  defined  by  CSX(NX),  CSY(NY), 
CSZ(NZ),  MULTEES  then  determines  the  8 grid  points  and  cor- 
responding grid  indexes  that  define  the  cell  in  which  the 
integration  point  (X,Y,Z)  is  located.  As  discussed  above, 
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however,  this  ceil  is  not  the  correct  prism  to  he  used  when 
interpolating  the  potentials.  Since  t he  curl  and  divergence 
of  the  displacement  field,  and  hence  the  potentials  them- 
selves, are  available  in  the  source  calculation  only  at  the 
centers  of  the  rectangular  cells,  MULTEES  must  next  deter- 
mine in  which  eighth  subcell  the  point  (.;,Y,Z)  is  located. 
Each  of  the  eight  subcell  regions  will  be  associated  with  a 
different  rectangular  prism,  in  which  the  eight  nodes  are 
centers  of  the  rectangular  grid  cells.  A 2-D  illustration 
of  the  method  under  discussion  is  presented  in  Fig.  G.l. 

In  this  example,  MULTEES  determines  that  the  point  P(X,Y) 
is  located  within  cell  (3,3),  that  is,  the  rectangular  cell 
bounded  by  the  grid  points  (3,3),  (2,3),  (2,2),  and  (3,2). 
MULTEES  then  determines  that  P is  located  in  the  "upper- 
right  subcell  region  of  cell  (3,5)  and  therefore  must  be  in 
the  rectangular  prism  formed  by  the  centers  of  the  four 
cells  (4,4),  (3,4),  (3,3)  and  (4,3),  where  a cell  (center) 
is  identified  by  its  "upper-righthand"  grid  point. 

Once  MULTEES  has  determined  the  rectangular  prism 
formed  by  the  eight  cell  centers  surrounding  the  point 
(X , Y , Z) , the  rectangular  coordinates  of  these  center  nodes 
are  computed  and  saved  along  with  the  cell  indexes  identify- 
ing each  of  the  eight  nodes.  In  this  manner  the  monitoring 
station  for  the  integration  point  (X,Y,Z)  is  determined. 

The  station  is  defined  by  the  three  indexes  corresponding 
to  the  grid  point  at  the  "origin"  or  "center"  of  the  eight 
cells  required  for  application  of  the  3-D  interpolation 
formula.  The  indexes  of  the  monitoring  station  are  identical 
to  the  cell  indexes  corresponding  to  the  seventh  node, 
following  the  labeling  convention  described  above.  For  the 
2-D  example,  illustrated  in  Fig.  G.l,  the  monitoring  station 
for  the  point  "P"  is  (3,3).  Notice  that  the  (3,3)  grid 
point  is  at  the  "center"  of  the  four  cells  whose  centers  are 
the  nodes  of  the  prism  in  which  P is  located. 
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Tlio  nr  XT  opcr.'i  i i on  pril'iinucd  bv  MIII.TM  i •.  i Ik  , I 
01  1 ;lt  ’ 0,1  ol  ,ll('  *'  * K*>  t ‘•liapc  ftiiR-l  i ons  corresponding  to  t In* 
cc" ter  nodes.  The  origin  of  the  normalized  coordinate  sys- 
tem is  given  by  Eq.  (G.28).  The  coordinate  transformation 
for  the  point  (X.Y,Z)  is  made  according  to  Eq.  (G.27).  Then 
the  shape  functions  S(p)  defined  in  Eq.  (G.25)  are  com- 
puted using  these  coordinates  (£,n,0  for  the  integration 
point  and  the  definition  of  the  eight  nodes  in  the  normalized 
coordinate  system  given  in  Eq.  (G.26). 


The  series  of  pre-processing  operations  described 
above  must  be  performed  for  each  integration  point 
In  order  to  calculate  the  appropriate  shape  functions,  a 
different  normalized  coordinate  system  is  computed  for  each 
point,  unless  possibly  more  than  one  point  is  located  within 
the  same  prism.  The  primary  output  generated  by  these  pre- 
processing operations  is  a set  of  eight  shape  functions  and 
a set  of  three  indexes  that  identifies  the  monitoring  station 
corresponding  to  each  integration  point.  The  acutal  appli- 
cation of  the  3-D  interpolation  formula  given  by  Eq.  (G.24) 
is  performed  by  a subsequent  program  in  MULTEEC  during  pro- 
cessing of  the  output  generated  by  the  source  calculation 
at  the  prescribed  monitoring  stations.  For  each  station, 
the  source  code  must  save  values  of  the  four  potentials  at 
the  centers  of  each  of  the  eight  cells  surrounding  the 
station. 
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G-7  SI LliCTION  OF  MONITORING  STATIONS 

Once  the  monitoring  station  corresponding  to  each 
integration  point  has  been  determined  as  described  above, 
three  further  pre-processing  operations  must  be  performed 
by  MULTEES  to  insure  efficient  linking  with  the  source  code. 
First,  the  set  of  monitoring  stations  are  sorted  by  MULTEES 
in  such  a manner  that  the  stations  will  be  arranged  accord- 
ing to  the  requirements  of  the  specific  source  code  to  which 
it  will  be  linked.  This  step  is  necessary  because  each 
source  code  operating  with  a rectangular  grid  system  pro- 
ceeds in  a prescribed  rectilinear  sequence.  The  monitoring 
stations,  on  the  other  hand,  have  been  determined  according 
to  the  angular  requirements  of  a spherical  integration.  By 
means  of  this  sorting  operation  MULTEES  can  arrange  the 
stations  to  suit  the  rectangular  sequence  followed  by  the 
source  code. 

After  sorting,  MULTEES  then  eliminates  all  duplica- 
tions among  the  set  of  monitoring  stations.  Depending  on 
the  angutar  interval  between  successive  integration  points 
and  the  linear  dimensions  of  the  rectangular  cells,  several 
integration  points,  especially  those  corresponding  to  small 
values  of  0,  may  be  located  within  the  same  prism  and 
therefore  be  associated  with  the  same  monitoring  station. 
While  the  shape  functions  for  such  points  will  be  different, 
only  a single  set  of  values  for  the  potentials  at  the  eight 
cell  centers  surrounding  the  monitoring  station  need  be 
sa^ed.  It  is  this  sorted  set  of  unique  monitoring  stations 
that  is  generated  by  MULTEES  and  passed  to  the  source  code. 

In  this  way,  the  source  code  is  informed  of  precisely  which 
cells  are  to  be  monitored  for  further  processing  by  MULTEES. 

As  the  third  step  in  the  selection  of  a set  of 
monitoring  stations,  a correspondence  table  must  be  pro- 
duced linking  each  integration  point  (I,J)  to  a sequence 
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number  lor  the  sorted  stations.  When  the  potentials  are 
passed  to  MUl.JT.hS  for  further  processing,  it  must  he  pos- 
sible to  determine  which  station,  identified  by  its  sequence 
number,  corresponds  to  any  given  integration  point.  The 
one-to-one  correspondence  initially  holding  between  a point 
(I,J)  and  a set  of  three  indexes  (IX, IY, IZ)  identifying  a 
station  must  be  preserved  even  after  the  stations  are  sorted 
and  the  duplications  are  eliminated.  This  association  is 
preserved  by  means  of  a correspondence  table  that  is  gen- 
erated and  saved  by  MULTEES  during  the  first  stage  or  pre- 
processing operation. 


PRELIMINARY  VALIDATION  OF  THE  MULTEES  PROGRAMS 


Preliminary  validation  of  the  MULTEES  program  series 
was  achieved  by  means  of  a comparison  with  the  multipole  co- 
efficients in  the  time  domain  calculated  previously  for  an 
underground  explosion  in  an  axisymmetric  cavity  and  reported 
by  Cherry,  e^t  al . (1974),  This  control  calculation  was  per- 
formed using  a two-dimensional  Lagrangian  finite  difference 
source  code.  The  divergence  and  curl  of  the  displacement 
field  required  for  the  equivalent  elastic  source  determina- 
tion were  monitored  at  twenty  evenly  spaced  points  on  an  arc 
subtending  90  degrees  at  a radius  of  140.68  meters. 


Two  approac *es  were  followed  for  testing  purposes.  In 
the  first  test  series,  by  applying  the  symmetry  properties  for 
an  axisymmetric  source  field,  a set  of  three-dimensional 


potentials  Xa(R»9><f>)  defined  at  evenly  spaced  points  on  the 


surface  of  a sphere  of  radius  140.68  meters  were  generated 
for  each  time  cycle  at  which  the  divergence  and  curl  were 
available  at  the  20  points  on  the  arc.  The  general  purpose 
MULTEES  program  was  then  used  to  compute  the  series  of  multi- 
pole coefficients.  Excellent  agreement  with  the  earlier 
results  of  better  than  1 percent  was  observed. 
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In  u second  test  series,  the  limited  set  ui'  values 
ol  xu  available  at  20  points  on  the  arc  corresponding  to 
a range  in  b of  0°  to  90°,  that  is,  output  from  a two- 
dimensional  source  calculation,  were  processed  directly  by 
MULTEES  under  the  2-D  mode  of  operation.  In  this  case, 
again  excellent  agreement  was  found  with  the  results  cal- 
culated previously. 

The  routines  for  applying  the  3-D  interpolation 
method  were  tested  separately  in  a number  of  ways.  In  one 
approach,  a set  of  analytic  functions  were  used  in  the  place 
of  values  of  the  potentials  at  the  eight  nodes  surrounding 
each  integration  point.  The  interpolated  value  at  each 
point  was  then  compared  to  the  exact  value  given  directly 
by  the  analytic  function.  In  this  manner,  the  accuracy  of 
the  interpolation  method  was  determined  as  a function  of  the 
position  of  the  integration  point  on  the  surface  of  the 
sphere  and  the  dimensions  of  the  rectangular  grid  system 
chosen  for  the  source  calculation. 
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